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Nonequilibrium mode-coupling theory (MCT) for uniformly sheared underdamped systems is de- 
veloped, starting from the microscopic thermostatted SLLOD equation, and the corresponding Liou- 
ville equation. Special attention is paid to the translational invariance in the sheared frame, which 
requires an appropriate definition of the transient time-correlators. The derived MCT equation 
satisfies the alignment of the wavevectors, and is manifestly translationally invariant. Isothermal 
condition is implemented by the introduction of the current fluctuation in the dissipative coupling 
to the thermostat. This current fluctation grows in the a-relaxation regime, which generates a devi- 
ation of the yield stress in the glassy phase from the overdamped case. Response to a perturbation 
of the shear rate demonstrates an inertia effect which is not observed in the overdamped case. Our 
theory turns out to be a non-trivial extension of the theory by Fuchs and Cates (J. Rheol. 53(4), 
2009) to underdamped systems. Since our starting point is identical to that by Chong and Kim 
(Phys. Rev. E 79, 2009), the contradictions between Fuchs-Cates and Chong-Kim are resolved. 

PACS numbers: 64.70.P-, 61.20.Lc, 83. 50. Ax, 83.60.Fg 



I. INTRODUCTION 

Liquids under shear are attracting continuous inter- 
est, not only for the importance of their application to 
industry, but also for their significance to the study of 
nonequilibrium statistical physics as ideal but experi- 
mentally accessible systems. Among them, special at- 
tention has been paid to dense supercooled liquids in the 
vicinity of the glass transition, i.e. glassy materials, due 
to their difficulty of understanding and their peculiarity 
compared with conventional systems which exhibit ther- 
modynamic phase transitions. 

Although not yet perfect, the mode-coupling theory 
(MCT) presents remarkable success in its application 
to glassy materials such as colloidal suspensions [1-3]. 
Probably the most striking feature of MCT is that it pre- 
dicts a two-step relaxation phenomenon (the /3-relaxation 
followed by the a-relaxation) [4, 5] of the intermediate 
scattering function (i.e. density time-correlator) in the 
vicinity of a critical packing fraction (p c , which is referred 
to as the "MCT transition point" . This two-step relax- 
ation is thoroughly investigated, and it is established that 
MCT is able to explain, e.g. the following properties: 
(i) square-root cusp anomaly of the temperature depen- 
dence of the plateau height of the density time-correlator 
(called the non-ergodic parameter, NEP) [6], (ii) power- 
law time-dependence with von Schweidler exponents of 
the density time-correlator around the /3-relaxation [7- 
9] , (iii) time-temperature superposition principle realized 
by the Kohlrausch- Williams- Watts-type behavior of the 
density time-correlator at the a-relaxation [10]. On the 



other hand, MCT is marred with the problem that the 
NEP can survive for temperatures below the critical tem- 
perature T c , while it decays to zero in actual glassy states 
[4, 11]. Moreover, ip c is about 10% smaller than the 
experimental value. To overcome these difficulties, suc- 
cessive studies beyond the conventional MCT have been 
carried out eagerly as well, e.g. (i) field-theoretic formu- 
lations [12-16], (ii) inclusion of higher-order correlations 
[17-19], (iii) estimation of the dynamic correlation length 
[20]. 

The framework of MCT itself is generic, and extensions 
to other materials such as dense granular materials [21- 
23] have also been attempted. In particular, MCT is 
able to incorporate the effect of shear and its resulting 
demolition of the "cage effect" . 

The introduction of shear into MCT has been worked 
out in the pioneering papers by Fuchs and Cates (FC) 
[24] and Miyazaki and Rcichman [25] , both of which have 
been formulated for sheared systems immersed in sol- 
vents, such as colloidal suspensions. The weak point of 
the approach of Ref. [25] , which was followed by Ref. [21] , 
has been pointed out by Chong and Kim (CK) [26]. 
They argued that the approach of Refs. [21, 25], where 
the steady-state structure factor (or, equivalently, the 
radial distribution function) is plugged in as an input, 
and then the steady-state shear stress is calculated as 
an output, is inconsistent, since they should be treated 
on the same footing. Rather, the equilibrium structure 
factor should be plugged in as an input, which is the 
case for Ref. [24]. (This scheme [27] is referred to as the 
"integration-through-transient (ITT)" in Ref. [28].) 

In this spirit, CK [26] have constructed MCT for a 
sheared system associated with a thermostat, where the 
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SLLOD equation [27, 29] is chosen as a microscopic start- 
ing point. In contrast to the theory of FC [24, 28, 30], 
where a Brownian system governed by the Smoluchowski 
operator is the starting point, the SLLOD equation is 
governed by the Liouville operator (Liouvillian), and the 
momentum variables are left unintegratcd, i.e. it is an 
underdamped system. 

Besides this issue, apparent discrepancies between the 
resulting equations of CK [26] and FC [24, 28, 30] have 
been noticed. For instance, the "initial decay rate" of the 
Debye relaxation, which exhibits the conventional Taylor 
dispersion, is time-dependent in FC [24, 28, 30], while 
it is not in CK [26]. Probably the most notable differ- 
ence is that there appears an additional memory kernel in 
CK [26], in addition to the conventional one which also 
appears commonly in other MCTs. As is well known, 
and is also argued in CK [26], these discrepancies can- 
not be the result of the difference of underdamped and 
overdamped systems, which exhibit quite similar behav- 
ior for long-time dynamics [31-33]. These contradictions 
were apparent, at least at the formal level, but as far as 
we know, no reconciliation has been proposed so far. 

In this paper, we reformulate the MCT for uniformly 
sheared underdamped systems, and resolve the contradic- 
tions mentioned above. We point out that the definition 
of the time-correlators by CK [26] should be modified to 
satisfy the translational invariance in the sheared frame 
and fulfill the requirement of physical correlations be- 
tween the Fourier modes of the fluctuations. This mod- 
ification leads to complexities, which are absent in the 
formulation of CK [26], due to the non-commutativity 
and the non-Hermiticity of the Liouvillians. The path 
which handles these complexities closely follows that of 
FC [30], which we believe to be the most sensible so far 
[34]. 

The crucial difference of our formulation to the CK's is 
the introduction of the isothermal condition, which is im- 
plemented at the level of the Mori- type equations. This 
is realized by incorporating current fluctuations into the 
dissipative coupling to the thermostat, accompanied by 
the introduction of a multiplier. The growth of this fluc- 
tuation in the a-relaxation regime generates deviation 
of the yield stress in the glassy phase from that of the 
overdamped case. The resulting equation formally corre- 
sponds to that of FC [30] in the overdamped limit, but 
the effect of the inertia or the fluctuating coupling to the 
thermostat renders this correspondence to be highly non- 
trivial. Hence, our framework appears to be a non-trivial 
extension of FC [30] to underdamped systems. The sig- 
nificance of the underdamped formulation is also demon- 
strated by its application to the calculation of the re- 
sponse to a perturbation of the shear rate, where the 
inertia effect is clearly observed. 

The paper is organized as follows. In section II, we 
briefly review the microscopic starting points. In section 
III, special attention is paid to the translational invari- 
ance in the sheared frame and the physically sensible defi- 
nition of the time-correlators. Adjoint Liouvillians are in- 



troduced, which turns out to be natural in the treatment 
of sheared systems. In section IV, Mori-type equations 
are derived by the application of the projection operator 
formalism, first without the isothermal condition. Af- 
ter that, the isothermal condition is formulated, and the 
Mori-type equation where this condition is implemented 
is derived. In section V, the mode-coupling approxima- 
tion is worked out, and a closed equation for the time- 
correlators (MCT equation) is derived. In section VI, a 
general formula of time-correlators for the steady-state 
quantities is presented, and the specific forms of the the 
dissipative coupling and the multiplier are discussed. A 
formula for the steady-state shear stress in this formula- 
tion is derived. In section VII, the result of the numerical 
calculation is shown for the density time-correlator and 
the steady-state shear stress, where the current fluctua- 
tion in the coupling to the thermostat plays an important 
role in the a-relaxation. The results of the CK theory [26] 
and the overdamped limit are also shown and are com- 
pared to the above result. In addition, response to a 
perturbation of the shear rate is shown. Section VIII 
is provided for discussion. Here, we first compare our 
work with the major preceding works. We hope that is- 
sues which were previously confusing are clarified. Then, 
the physical significance of the present formulation and 
the possibility of its extension are discussed. Finally, 
in section IX, we summarize our results and conclude 
the paper. An application to the response theory (Ap- 
pendix A) , comparison of the Mori-type equations in the 
underdamped and overdamped cases (Appendix B), and 
technical details (Appendices C and D) are collected in 
the Appendices. 



II. MICROSCOPIC STARTING POINTS 

Here we briefly summarize our microscopic starting 
points, i.e. the SLLOD equation, the Liouville equation, 
and the steady-state formula. Our treatment closely fol- 
lows that of CK [26] , so the details in common arc to be 
referred to Ref. [26]. The crucial difference, however, is 
that our formulation is intended to satisfy the isother- 
mal condition, which will be presented in section IV D. 
Attentions related to this issue will be paid. 



A. SLLOD equation 

We deal with an assembly of N equivalent spheres with 
diameter d and mass to, interacting with themselves and 
a thermostat in a volume V . The interaction between the 
spheres is assumed to be a two-body soft-core potential 
force. Uniform shear is imposed on this system, where 
the shear velocity is given by = n ■ r. Here, k is the 
time-independent shear-rate tensor, which is assumed to 
be of the form k^ v = jS^S^y in this paper, where the 
shear rate is denoted as 7. r is the spatial coordinate, 
and Greek indices \,(j,,v,--- denote spatial components 
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{x, y, z} in the remainder. We introduce the distance of 
the two shear boundaries, L, for later convenience. Then, 

the shear velocity at the boundaries are = ±jL/2 
for y = ±L/2, respectively. The Newtonian equation of 
motion for the i-th sphere (i = 1, 2, • • • ,N) is given by 
the following SLLOD equation [27, 29], 



Pi(t) 



pM 

m 



+ K-n(t), (l) 

-K-pi(t)-a(r) Pi (t). (2) 



Here, , Pi{t)} is the position and the momentum of 

the i-th sphere at time t, and T(t) = , Pi{t)}jL\ 

is the phase-space coordinate. In Eq. (1), Pi(t) is de- 
fined as a relative momentum in the sheared frame, and 
is referred to as the peculiar or the thermal momentum. 
Fi(t) = -dU/dn(t) = -£,• 0u ('«(*)) /0r«(i) is the 
conservative force acting on the i-th. sphere from other 
spheres, U is the total potential, u(rij(t)) is the two- 
body potential, (t) = n (t) — rj (t) is the relative posi- 
tion and rij(t) = \rij(t)\ is the relative distance, between 
the i-th and j-th spheres, respectively. The parameter 
a (r) represents the strength of the coupling of spheres 
to the thermostat, which in general implicitly depends on 
time through its dependence on T(t). Note that ot(T) has 
been assumed to be a constant in Ref. [26] . It is possible 
to prevent the system from heating-up and control the ki- 
netic temperature to be a constant by tuning the coupling 
a(r), in which case a(T) should be regarded as a multi- 
plier, rather than an independent physical parameter. A 
well-known example can be found in molecular dynamics, 
which is the Gaussian isokinetic (GIK) thermostat [27]. 
This thermostat is distinguished by its generality and its 
importance, which has the explicit time dependence as 



a(t) 



E t [F i (t)-p i (t)--ytf(t)pV(t)} 



(3) 



which can be derived from the condition 4? £jPi(£) 2 = 
0. The specific form of a(T) in MCT with isother- 
mal condition will be discussed later in section VI B. 
Note that Eqs. (1) and (2) reduce to Newtonian equa- 
tions if we introduce the viscous force Fj vis \t) = 
—a(T)m (ri(i) — k • ri(i)) , except for the instant onset 
of the shear. 



B. Liouville equation 

The Liouville equation for a phase-space variable 
A(T(t)) reads 



r=r(t) 



e lL HLA (T(0)) = i£A (T(t)) , (4) 



where the action of the Liouvillian iC is defined as 

dA{V) 



iCA(T{0))^T{T) 



dV 



(5) 



Here, and in the remainder, the abbreviated notation 
r = r(0) is adopted. The formal solution of Eq. (4) is 



given by A (T(t)) 



jet 



A (T). Note that, in Eqs. (4) and 



(5), the Liouvillian is assumed not to bear explicit time- 
dependence. This is compatible with a time-independent 
shear, which we consider in this paper. 

On the other hand, the Liouville equation for the 
nonequilibrium distribution function p (T, t) reads 



d P (T,t) 
dt 



p(T,t) = -itfp(T,t), (6) 



where 



i£t = f • — + A(r) = ic + A(r) 



is the adjoint Liouvillian, and 



A(T)^-f 



(7) 



(8) 



is the volume contraction factor of the phase space. The 
formal solution of Eq. (6) is p(T,t) = e~ ich p{T,Q). In 
general, A(T) ^ for nonequilibrium systems, and hence 
the Liouvillian is non-Hermitian. In our case, A(T) is 



= -3JVc,(r)-yj „ 



Opt 



-3Na(T), 



(9) 



where the last approximate equality holds in the ther- 
modynamic limit. This will be verified for our specific 
choice of a(T) in section VI B. 

Let us decompose the Liouvillian as follows for conve- 
nience: 



iC = iCq + iCy + iC a , 

?i • — + F t ■ — 

m dri 1 dpi 



ic = Y J \- — i' ■ 



(10) 

(11) 



iC a = -a(T)^2pi 



d 

dpi 



(13) 



r'=r 



As can be easily seen, iCo, iCj, and iC a are the time- 
evolution generators of the conservative force, shearing, 
and the thermostat, respectively. We consider a situation 
where the system is initially in equilibrium with temper- 
ature T, and at time t = a uniform shear with shear 
rate 7 and a thermostat with strength a(T) is turned on. 

The Liouvillian obeys the following adjoint relation in- 
side the integral with respect to the phase-space coordi- 
nates, 

J dT [iCA{T)] B{T) = - j dVA{T) [iC^B{T)] , (14) 
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whose repeated use leads to 
J dT [e lCt A(T)] B(T) = J dTA(T) 



B(T) .(15) 



Equation (14) can be easily shown by partial integration. 
The adjoint relation Eq. (14) and the definition of the ad- 
joint Liouvillian Eq. (7) leads to the following important 
relation which indicates the action of the Liouvillian in- 
side the time- correlators, 

([iCA(T(t))]B(Tr) 
= -(A(T(t)) [iCB(T)]*) + (A(T(t))B(T)*n(T)) ,(16) 

where f2(r) is the work function, 

Q(T) ee -pia xy (T) - 20a(T)5K(T). (17) 

Here, a xy {T) = £i (pfPi /m + Vi F?) and 5K(T) = 
K(T) — ZNksT/2 are the zero-wavevector limit of the 
shear stress and the fluctuation of the kinetic energy 
^( r ) = E 4 P?/( 2m ), respectively, and /3 = 1/ (k B T). 
The derivation of Eq. (17) is shown in Appendix CI. 
The ensemble average (• • •) in Eq. (16) is defined for a 
phase-space variable A(T(t)) as 

(A(T(t))) = J dT Pini (T)A(T(t)), (18) 

where p[ n i(T) is the initial Maxwcll-Boltzmann distribu- 
tion, and the definition of Eq. (18) corresponds to the 
"Hcisenbcrg picture" [27], which we adopt in this work. 

Finally, we state here a steady-state formula for a 
phase-space variable A(T(t)), 



poo 

A) SS = (A(T))-Pj / dt(A(T(t))a xy (T)) 
Jo 



-2/3/ dt(A(T(t))a(T)5K(T)}, (19) 
Jo 



where 



(A) ss ee lhn (A(T(t))) 



(20) 



is the ensemble average in a nonequilibrium steady-state. 
Refer to Appendix C 1 for the derivation of Eq. (19). 



iq(-t)-f 



A. Fourier transform 

A crucial feature of the sheared system is that the 
translational invariance is preserved only in the sheared 
frame. Hence, in order to examine the implications of the 
translational invariance, we must move on to this frame. 
As derived in Appendix C2, Eq. (C22), the wavevector 
of the Fourier transform in the sheared frame is Affinc- 
deformed from that of the experimental frame (we will 
refer to it as the "advected wavevector" in the following): 

A q( _ t) {t) = j d 3 fA(f,t)t 

q(t) ee q — q ■ nt. 

Here, t = t and f = r — (k ■ r)t are the temporal 
and spatial coordinates in the sheared frame, respec- 
tively. We adopt the definition of the advected wavevec- 
tor of FC [30] , which differs from the conventional one, 
q(t) ee q + q-Kt [21, 25, 26]. The reason of this choice will 
be explained later. The values of the phase-space vari- 
ables are equivalent, irrespective of the frames, so the 
following equality holds: 

A q( _ t) (i) = A q {t). (23) 

The time-evolution of the Fourier transform in the 
sheared frame is generated by iCo, iC a , and i£^ p , where 
i£-y p is the momentum part of iC^: 



(21) 
(22) 



—-A q( _ t) (t) = iCA q( _ t) (t), 



iC ee iC + iC a + iC 



iC 



d 



ET u 



(24) 
(25) 
(26) 



The derivation is shown in Appendix C2. Note that the 
action of the coordinate part of i£^, which we denote 



K 



dri 



(27) 



is already incorporated in the advected wavevector. This 
fact already has been pointed out by Hayakawa and Ot- 
suki (HO) [21]. We refer to i£j r as the "advection Liou- 
villian" , since it generates an advection of the wavevec- 
tors of plane waves, e ~ lC ^ rt e lq r — e lq ^' r . 



III. TRANSLATIONAL INVARIANCE IN THE 
SHEARED FRAME 

So far we have repeated the formulation by CK [26], 
aside from the introduction of the dependence on the 
phase-space variables T of the dissipative coupling to the 
thermostat, a(T). In this section, a detailed discussion 
is devoted to the translational invariance in the sheared 
frame, and to the definition of the time-correlators. 
This discussion clearly indicates a crucial mistake in the 
Fourier transform of the phase-space variables in CK. 



B. Adjoint Liouvillians inside the time-correlators 

We have seen above that the Fourier transform sepa- 
rates the Liouvillian iC into iC and the advection Liou- 
villian iC^ r . Hence it is convenient to decompose the ad- 
joint relation Eq. (16) into those for iC and iCj r , and de- 
fine the corresponding adjoint Liouvillians, i£) and i£\ r - 
As for iC, it is 



iCA(T(t)) B(T) 



A(T(t)) \ifrB(T) 

iC — 0, 



,(28) 
(29) 
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where 



n(T) = -Pia^ (T) - 2Pa(T)5K(T) (30) 



is the modified work function, which includes only the 
kinetic part of the shear stress, 



(kin) 
'xy 



E 



m 



As for i£j r , it is 



([i^ r A(T(t))]B(ry)= (A(T(t)) \iC\MT) 



iC 



t - 



(31) 



,(32) 
(33) 



where the repeated use of Eq. (32) results in 
( [e"^r*A(r(t))] B(T)*) 
= (A(T(t)) 



e - i£ l*S(r) 



(34) 



Here, <r xy 



( - pot - ) is the potential part of the shear stress, 

^f^E^' ( 35 ) 



Note that the adjoint Liouvillians iC) and i£\ r are well- 
defined only inside the time- correlators, and are not to 
be confused with the adjoint Liouvillian iC\ which is 
defined in Eq. (7) as an independent operator. 



C. General time-correlators 

In nonequilibrium statistical mechanics, time- 
correlators play an essential role. Hence, we figure out 
the implications of the translational invariance on the 
time-correlators. An immediate consequence of the 
translational invariance in the sheared frame is 



A(f,t)=A(f + a(t),t), 
A k (t) = e ik - a ^A k (t), 



(36) 
(37) 



where a(t) = (jtL, L, 0). Here, L is the distance between 
the two shear boundaries, which is introduced in section 
II A. Equation (37) leads to the following "selection rule" 
of the wavevectors of the one-point and two-point func- 
tions in the sheared frame: 

(A k (t)) = {A k=0 (t))S kfl , (38) 
(A k (t)B* q (0)) = (A q (t)B* q (0))5 k , q . (39) 

Application of the equivalence of the Fourier trans- 
forms in the experimental and the sheared frames, i.e. 
Eq. (23), to Eq. (39) results in the following "selection 
rule" in the experimental frame: 

(A k{t) (t)B*(0)) = (A q{t) (t)B*(0))S k , q . (40) 



Here, the Fourier transform with an advected wavevector 
A q (t)(t) is explicitly written in terms of the Liouvillians 

as 

A q{t) (t) = ^ iCt A i (T(0))eW"' 

i 

= z lCt ^ ( r (°)) e- iC *<-*e iq - rt , (41) 

i 

where Ai (T(t)) is a Fourier coefficient which is defined 
in Appendix C2, Eq. (C9). 

One might think that the two-point function Eq. (40) 
involves an ambiguity; it might seem that another choice, 

e.g. ^ A q (t)B*^ _ ( )(0) ^, is equally valid. Actually, 

this has been the choice made in CK [26] and FC 
[24, 28]. However, the two expressions (A q ^(t)B*(0)) 

and ^ A q (t)B*^ _ t j(0) ^> arc in fact incquivalcnt, due to 

the non-commutativity and the non-Hcrmiticity of the 
Liouvillians. We will prove this statement below. First, 
with the use of Eq. (41), ( A q ^(i)B*(0) ) can be written 
in the following form: 

(A q{t) (t)B* q (0)) 



J2^ Ct A l (T(0))e- lC ^ t e t<1 - r ' 



iq-Vj 



(42) 



Similarly, ^ A q (t)B*^_ t ^(0) ) can be written as follows: 
A q (t)B* q{ _ t) (0) 



J2e lCt A t (r(0))e lq - r ' 



5>,(r(0))V £ - 



(43) 



Even when the "advection Liouvillian" i£j r commutes 
with Ai(T(0)) and Bj(T(0)), which is the case of interest 
in MCT where these variables are the density and the 
current-density fluctuations, Eqs. (42) and (43) are not 
equivalent. This can be seen by the use of the adjoint 
relation of the Liouvillians Eq. (34) and the relations 
[i£^ r ,i£] ^ and iC^ r ^ i£t . It can also be fore- 
seen from Eqs. (42) and (43) that different definitions of 
the two-point functions lead to different physical conse- 
quences. 

We assert that the specific definition Eq. (40) is the 
physically sensible choice. It states that a fluctuation at 
time t = with a wavevector q = q(0) is correlated at 
time t with a fluctuation with a wavevector q(t), exclu- 
sively. Note that the definition of the advected wavevec- 
tor Eq. (22) is chosen for the compatibility to the intuitive 
picture described above. This definition is essentially co- 
incident with the one adopted in the previous studies 
[21, 25, 30], although it appears to be ( A q{ _ t) (t)B*(0) ) 
in Refs. [21] and [25]. This superficial discrepancy with 
Eq. (40) is due to the different definition of the advected 
wavevector, q(t) = q + q ■ nt. 
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IV. MORI-TYPE EQUATIONS 

In liquid theory, slowly-varying (i.e. long-wavelength 
and low-frequency) conserved variables are of inter- 
est. Conventionally these are the density fluctuation 
and the current-density (momentum) fluctuation. In 
MCT for sheared thermostatted systems, the formula 
for the steady-state quantities is formulated in terms of 
a time-correlator for density fluctuations (density time- 
correlator). In this section, we introduce the time- 
correlators of interest and derive the Mori-type equations 
[35] for them by applying the projection operator for- 
malism. We first derive the Mori-type equations without 
isothermal condition. Then, we formulate the isothermal 
condition by introducing a multiplier for this constraint. 
Finally, we derive corrections to the Mori-type equations, 
where the multiplier is introduced. 



A. Time-correlators of interest 

Density and current-density fluctuations at equilib- 
rium are denoted in Fourier space as 



(44) 



Jq 



J2-e iqri (X = x,y,z), (45) 
jL — 1 m 



space. Note that U(t) is non-unitary since the Liouvil- 
lians are non-Hcrmitian. Then, Eqs. (46) and (47) are 
expressed respectively as 

*«(t) = ±([U(t)n q ]n' q ), (51) 

H fr) = Jj (52) 



B. Continuity equations 

Now we derive the equations of motion for Eqs. (51) 
and (52). The time derivative of U(t) is 



d_ 



U(t) = e lCt (iC-iC^)e 



= e lC HCe 



, (53) 



where iC is defined in Eq. (25). The action of iC on n q 
is obtained from Eqs. (11), (13), and (26) as 



i£n q = iq ■ j q . 
From Eqs. (52)— (54) we obtain 



(54) 



d . , 



d rT , \ 
-U(t)n q 



i iCt iCn 



<?(*) 



q(t)-H q (t). (55) 



respectively. The spatial dimension is assumed to be 
three, in accordance with the numerical calculation which 
will be carried out in section VII. The corresponding 
time-correlators of interest of the form Eq. (40) are as 
follows: 



%(t)^-(n q{t) (t)n q (0y) 



^(^-(4,(^,(0) 



(46) 
(47) 



Here, $ q (t) is the density time-correlator, and H^(t) is 
referred to as the density-current cross time-correlator. 
As already mentioned, the Fourier coefficient Ai(T(0)) 
in Eq. (41) is Ai(T(0)) = 1 for the density fluctuation 
and Ai(T(0)) = p^/m for the current-density fluctuation; 
hence they commute with iC^ r , respectively. This leads 
to the following expression for the density and current- 
density fluctuations at time t, 



£,(*)(*) = = e^% (t) = U(t)^ q , (48) 



where 



U{t)=e lCt e- lC ^ 1 . 



(49) 
(50) 



Here, £ is cither of the hydrodynamic variables, n or j, 
and U(t) is the time-evolution operator for £ in Fourier 



On the other hand, the action of iC on j q cannot be 
written in a concise form as simply as Eq. (54): 



i 

= N 



dt 



(56) 



A conventional way to handle Eq. (56) is to deform it 
into a Mori-type equation [35] . This task is conducted in 
the next subsection by the application of the projection 
operator formalism. 



C. Projection operator formalism 

We have introduced the time-correlators, Eqs. (46) and 
(47), in a physically sensible way that a fluctuation at 
time t = with a wavevector q is correlated at time t 
with a fluctuation with a wavevector q(t). We refer to 
this feature as the "alignment of the wavevectors" in this 
paper. Even if the wavevectors of the time-correlators 
are aligned, this feature is not necessarily preserved in 
their entire continuity equations. As for the density time- 
correlator $q(t), this is positive as can be seen from 
Eq. (55). In deriving a Mori-type equation for the cross 
time-correlator H^(t), we demand the alignment of the 
wavevectors as a principle. 
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For this purpose, we introduce the following time- 
dependent projection operator [30] 



Xn 



fe(t) 



NS, 



nk(t) 



E 



k - jfe « k 

and its complementary operator 

Q{t) = l-V(t). 



(58) 



Here, X is an arbitrary phase-space variable in Fourier 
space and the normalization factors are determined from 
the equal-time (equilibrium) correlators, ( n q n*, ) = 

NS q S q ^ q ' and (jqjg,^ = Nv%6 x » 6 q , q < , where S q is the 

static structure factor [35] and vt = \JksT jm is the 
thermal velocity. These operators preserve the desired 
properties, (i) idempotency: V(tf = V(t), Q(t) 2 = Q(t), 
(ii) orthogonality: V{t)Q(t) = Q(t)V(t) = 0, and (hi) 
Hcrmiticity: 

( [P(t)A q(t) {t)} B q (0y) = (A q(t) {t) [V(t)B q (0)r), (59) 
< [Q(t)A q{t) (t)] B q (0y ) = (A q{t) (t) [Q(t)B q (0)r ) . (60) 

In addition, we further introduce a "rescaled static pro- 
jection operator" [30], 



p t x = j: 



NS k 



W>,> 



(61) 



whose raison d'etre will be explained later. Although 
Eq. (61) is a projection operator onto the subspace 
spanned by the static density and current-density fluctu- 
ations {rife, jfc}, it involves a dependence on time through 
the advected index of S^t) , so we appended a subscript t. 
One can easily verify, by the_use of Eq. (34), the following 
relation between V(t) and P t , 



V(t) = e- lC ^V t e 



iC\ t 



(62) 



where the adjoint of the advection Liouvillian iC!y r is 
defined in Eq. (33). 

Now we derive a Mori-type equation for Eq. (56) by 
inserting the projection operators as follows: 



dt 



U(t) = e lCt [P(t) + Q(t)} iCe 



iCt 



iCe 



= e iCt e~ iC ^> rt 



V t e 



iCe 



U(t) 



where Uo(t,t') is the solution of the homogeneous 
equation. Uo(t,t') can be written in terms of the 



time-ordered exponential, exp_ 



J, dsX( S ) 



= 1 + 



ife(t), (57) E^=i /** ds i • • • It" 1 ds nX{s n ) ■ ■ ■ X( Sl ), as follows: 



(63) 



The formal solution of Eq. (63) is given as 

U(t) = U (t,0) + f dsU(s)P s e iC l s iCe- lC ^ s U {t,s), 
Jo 

(64) 



U (t,t') = exp_ 



dse lL ^ s Q{s)iCe 



(65) 



From Eqs. (63) and (64), the time derivative of the 
current-density fluctuation, which appears in Eq. (56), 
can be decomposed into the "correlated part" and the 
"uncorrelated part" : 

d 



U(t)j x = U{t)V t e lC ^iCj x q(t) + Uo(t,0)e iC ^Q {t)i l^ 



L 



+ / dsU(s)P s e l ^r s iCe~ lL ^ s U {t,s) ■ e 



(66) 



Here, the first term on the right-hand side (r.h.s.) is the 
"correlated part" and the second and the third terms are 
the "uncorrelated part". Applying Eq. (61) to Eq. (66), 
and then substituting Eq. (66) to Eq. (56), we obtain 
the Mori-type equation. Following the derivation in Ap- 
pendix C3, we obtain 

q(t) x 



1^' 



-v. 



S Q(t) 



N 



+ ~ ( Rq {t) (t)n q 



f 

Jo 



d S L X (t, S )<i>q( S ) 



dsM^(t,s)H%(s), 



where we have introduced 



^(t)(*) 

R q(t) 



iL x q {t,s) 



U Q (t,0)e lC ^ t R x 



<?(*)' 



i 



NS qW 
1 



i£U (t 7 s)R x 



■«(*) 



ijCU (t,s)R q{t) 



(67) 

(68) 
(69) 

(70) 
(71) 



U (t,s) = e- lC ^ s U (t,s)e 



AC-, 



(72) 

Here, Rq^it) is the "random force", whose time- 
evolution is given by the projected time-evolution op- 
erator, Uo(t,0). There appear two types of memory ker- 
nels, L x (t,s) and M^(t,s), due to the projection onto 
the current-density and the density fluctuations, respec- 
tively. The correction of the friction coefficient ao due 
to the introduction of the isothermal condition will be 
discussed in section IV E. 

Note that the "random force" Eq. (68) is not orthogo- 
nal to the density fluctuation at this stage. An additional 
requirement leads to the orthogonality, which will be dis- 
cussed in Eqs. (88) and (89) in section V. Note also that 
the memory kernels possess two time arguments, in con- 
trast to those which appear in CK. This issue will also 
be discussed in section V. 
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D. Isothermal condition 

So far we have derived the Mori-type equations with- 
out isothermal condition, which is a condition necessary 
to hold the kinetic temperture unchanged. In this subsec- 
tion, we attempt to implement the isothermal condition. 
The derivation of the resulting Mori-type equations will 
be performed in section IV E. 

At the level of the SLLOD equation, it is known that 
the isothermal condition is satisfied by the specific choice 
of a(t) given by Eq. (3), which is referred to as the GIK 
thermostat. In principle, the Mori-type equations de- 
rived from the corresponding Liouville equation, with 
a(r) given by Eq. (3), no longer contain any multiplier, 
and the invariance of the kinetic temperature is assured 
automatically. However, this choice of a(T) leads to a 
difficulty in deriving the Mori-type equations, since the 
integral of a rational function with Gaussian weight is 
difficult to perform explicitly. 

This motivates us to implement the isothermal con- 
dition at the level of the Mori-type equations, which is 
attained by requiring the time derivative of the average 
kinetic temperature to vanish. Note that this scheme re- 
quires the time derivative of the kinetic temperature, not 
necessarily itself but at least its average, to vanish. In 
this scheme, we avoid to fully incorporate the fluctuations 
as in Eq. (3), and attempt to incorporate them partially, 
compensating it with a multiplier at the level of the Mori- 
type equations. We denote this multiplier, which corre- 
sponds to the multiplier of the Liouville equation a(T), 
as X a (t) in the remainder. 

By differentiating the generalized Green-Kubo formula 
for the kinetic temperature, we obtain 



!<T(r ( t))> 



3Nk B 

2/3 



3Nk B 



(K(T(t)MT)) 
{jG K!a (t) + 2G K ,aSK(t)}, (73) 



where 



and 



G K , a (t) = ([U(t)K(T)}a xy (T)) 



(74) 



a(r). This leads us to define a "renormalized" time- 
evolution operator U R (t), where 



U R {t) = X a {t)U{t), 



(78) 



and define G^ aSK (t) as 



G'Z SK (t) = ([U R (t)K(T)]a(T)SK(r)) 

= X a (t)GK,aSK(t)- 



(79) 



Here, X a (t) factors out, since it is not a phase-space vari- 
able. Then, if Gic,a8K(t) ^ 0, choosing the multiplier 

as 



A a (t) 



7 GkA*) 

2 GK,aSK(t) 



(80) 



assures the isothermal condition Eq. (77) to be satisfied. 
Note that, while a(T) is expressed in terms of the phase- 
space variables T, X a (t) is expressed in terms of time- 
correlators at time t, when the constraint Eq. (77) is 
satisfied. Concrete examples of a(T) and X a (t) will be 
discussed in section VI B. 



E. Isothermal Mori-type equations 

Now we figure out how the multiplier X a (t) enters in 
the Mori-type equations. The identification of the com- 
plete correction is possible, but rather lengthy. Further- 
more, most of the modifications vanish in the mode- 
coupling approximation, which we will discuss in sec- 
tion V. Hence, we only focus on the modification which 
survives the mode-coupling approximation, which is for 
the friction term a H^(t) in Eq. (67). As shown in Ap- 
pendix C 3, this terms arises from the "correlated part" 
in Eq. (66), which is explicitly 



U(t)P t e 



<?(*) 



U{t)P t e iC ^a{T)iC v j^ t) 



(81) 



Here, we extract a(T) and define the remaining part as 



G K , aSK {t) = < [U(t)K(T)}a(T)6K(T) ) . (75) done in Eq. (78), which leads to 



J2iPi ' 9/dpi. We introduce \ a (t) here as is 



Hence, the isothermal condition at the level of the Mori- 
type equations reads 

iG Kta (t)+2GK,aSK(t)=0. (76) 

In order to satisfy this constraint, we introduce X a (t) in 
GK,aSK{t)i the result of which we denote G^ aSK (t), and 
recast Eq. (76) as 

jG K At)+2G% aSK (t)=0. (77) 

However, from Eq. (75), we can see that the definition 
of G^ aSK (t) is non-trivial, since the effect of the time- 
evolution of the multiplier appears in U(t), rather than 



<) = X a {t)a 0l (82) 



U R (t)V t e iC ^a(r)iC p j^ t) 



where the explicit form of a(T) which will be specified in 
section VI B, Eq. (112), is utilized. Hence, the Mori-type 
equation with isothermal condition is 



d_ 

Jt~ q 



H X q (t) = -v 2 T ^$ q (t) - X a (t)a Q H$(t) - [k ■ H q (t)Y 



(83) 



+-( Rq - {t) (t)n* q )- dsL$(t,8)* q (8) 



f dsM^{t,s)H»{s). 
Jo 
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The fact that X a (t) does not appear in the memory ker- 
nels in the mode-coupling approximation can be verified 
in section V, where the resulting memory kernels turn 
out to be independent of a®. In addition, as shown in 

section V, the vanishing of ^ R X ^ (t) n q ^ is not affected 
by a(r), and hence X a {t) does not appear here either. 

In principle, the Mori-type equations which follow from 
the choice Eq. (3) of a(T) contain no multiplier, and is 
expected to satisfy the isothermal condition Eq. (7G) au- 
tomatically. On the other hand, the Mori-type equations 
Eqs. (55) and (83), which follow from a(T) of Eq. (112), 
contain a multiplier X a (t), and the isothermal condi- 
tion Eq. (77) is explicitly satisfied if X a (t) is chosen as 
Eq. (120), which will be derived in section VI B. We ex- 
pect that the physical equivalence is achieved for these 
two schemes. 

The Mori-type equations of the underdamped and 
overdamped [30] cases are compared and discussed in Ap- 
pendix B. 



MODE-COUPLING APPROXIMATION 



and (71), iCUo(t,s)R x ^, can be deformed as follows, 



i£U (t,s)R x 



•«(*) 

it [1 + £(*)] e- lC l s U (t, s)e"^*<2( t )flA (t)) (85) 



where Eq. (72), the idempotency of Q(t), and the identity 



[1 + £(*)] e" 



■iCl t 



is applied. Here, 

£(i) = 07 /' 
Jo 



dse 



iCt s (pot) „-iC^ 



a xy ' e 



(86) 



(87) 



is the accumulated clastic energy due to shear, where 
cTxy^ is the potential part of the shear stress introduced 
in Eq. (35). As discussed in FC [30], the shear- induced 
term £(i) is an obstacle for the application of the second 
projection operator and the factorization approximation. 
We assume here £(i) ~ 0, whose validation is discussed 
in FC [30]. At least, this assumption is valid in the weak- 
shear regime, since £(£) is proportional to 7. 

The neglect of £(£) leads to the following relation, 



We have derived an isothermal Mori-type equation for 
the cross time-correlator H x (t), Eq. (83), in the previ- 
ous section. However, this is not a closed equation for 
the time-correlators $ q (f) and H x (t), unless the memory 
kernels are expressed in terms of them. For this purpose, 
we introduce a time- dependent second projection opera- 
tor [30] , which extracts the dynamics correlated with the 
slowly-varying pair-density modes: 



V 2 (t)X EE 



k> P N2S Ht) S p(t) 



nk(t) n P (t)- 



(84) 



The normalization factor is determined by the fac- 
torization approximation of the equal-time (equilib- 
rium) four-point function of the density fluctuations, 



Sk,k>Sp,p>N 2 S k ( t )S p ( t ) (k > p,k' > p'). The second 
projection operator Eq. (84) is idempotent and Hcrmi- 
tian, similar to the projection operator Eq. (57). 

In underdamped MCT, there exists a potential next- 
to-leading second projection operator, which projects the 
dynamics onto the density-current modes [36]. How- 
ever, it can be shown that this projection is negligi- 
ble as compared to the projection onto the pair-density 
modes, at least for time-reversible thermostatted systems 
[37]. Hence, the choice of Eq. (84) is assured. Note 
that this feature does not always hold, e.g. for granular 
systems, where the projection onto the density-current 
modes plays an essential role [23, 36]. 

There is one subtle issue we should handle in order for 
the application of Eq. (84) to work [38]. The operator 
which appears in the memory kernels defined in Eqs. (70) 



e- lC l s U (t,s) = Q(s)e- tC l s U (t, S ), 



whose proof is given in Appendix C 4. The first implica- 
tion of Eq. (88) is the orthogonality of the random force, 



U Q {t 1 {))e lC ^ t R x 



Q(0)U o (t,0)e iC ^ t R x 



q(t) 



U o (t,0)e lC ^ t R x 



i(t) 



Q(0)C ) - 0, (89) 



where £ = n or j. The second implication is the form of 
the aforementioned operator described in Eq. (85), 



iCU (t,s)R x {t) 



~ iCQ(s)e-^U (t,s)e^tQ {t ) R ^ ( 90 ) 

which now has the desirable feature, i.e. the projection of 
Uo(t,s) is now complete. The memory kernels Eqs. (70) 
and (71) can be rewritten, from Eq. (90), as follows: 



iL x (t,s) = 



1 



NS, 



<?(*) 



m, S )R x 



Q(s) 



nq(s)^ 



M x »(t, ; 



U' {t,s)R x 



<?(*) 



9(a) 



Here. 



U^(t,s) ee Q(s) e - iC l s U (t,s)e iC ^ t Q(t) 
is the modified projected time-evolution operator, 



3q{s) U 



,(91) 



(92) 



(93) 



(94) 
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is the modified random force, and is the modified work 
function defined in Eq. (30). The derivation of the above 
equations Eqs. (91)-(94) are shown in Appendix C 5. 

Now we insert the second projection operator Eq. (84) 



into Eqs. (91) and (92) as Q(s)e~ tC l- s Uo(t, s)e iC ^ t Q{t) 

V2(s)Q(s)e- lC ^ s Uo{t,s)e lC ^ 4 Q(t)P 2 (t), which re- 
sults in the following form, i.e. products of vertex func- 
tions at times s and t, bridged by a propagator from time 
s to t: 



iL x q {t,s) 



NS, 



9(t) 



EE- 

k'>p' k>p 



R q(t) n *k(t) n *p(t) 

N 2 Sk(t)S P (t) 



EE 



Nvl — — 

1 k'>p' k>p 



R q(t) n k(t) n *p(t) 

N2 Sk{t)Sp(t) 



Uo(t, s)n k{t) n pW 



n k'(s) n p'(s) 



Uo(t, s)n k ( t )n p (t) 



n 



k'( S ) n p'{s) 



N 2 Sk>( s )Sp'(s) 



n k'(s)n p '( s )AR^ s) 



(95) 



(96) 



We can see from the above expressions that the derived 
forms of the memory kernels are consistent with the prin- 
ciple of the "alignment of the wavevectors" ; the vertex 
function at time t includes as indices only the advected 
wavevectors with argument t, e.g. q(t), and a similar 
feature also holds for the propagator. 

The remaining tasks are the calculation of the vertex 
functions and the approximation of the propagators . As 
for the vertex functions, the convolution approximation 
[35] is applied. We only show the results below, since the 
derivation, which is shown in Appendix C6, is straight- 
forward. 



DA * * 

n q(t) n k(t) n p(t) 

N 2 S k ( t )S p (t) 



n k{t) n p{t) AR q * (t) 



q,k+pVq(t),k(t),p(t) 



(97) 



■ n x 



-jj°q,k+ P V£ t)Mt)Mt) , (98) 
(99) 



A* 



N 2 S k (i)S p (t) 

V q X k ,p = v 2 (fc A c fe +Ap) 



n k (t)n p ( t )Q(t) 



n q(t) n 



= 0. 



(100) 



As for the propagator, we adopt the factorization ap- 
proximation, which replaces it with the product of the 
projection-free propagators: 

U' a {t,s)n k{t) n p(t ) n* k , (s) n pl {s) 

- <W,feV,P $ fc( S )(* - »)* P (a)(* - s). (101) 
The derivation of Eq. (101) is shown in Appendix C 7. 

From Eqs. (95)~(101), we arrive at the final expressions 
for the memory kernels, 

iL q (t,s)=0, (102) 



A 



T J (2tt) 

X$k(s)(t - S)$ p{s) {t 



i ■ q(t)Mt),P(t) q(s) Ms) Ms) 

S), (103) 



where the summation of the wavevectors is replaced by 
the integral, and p = q — k. Note that the memory ker- 
nel L q introduced by CK [26] vanishes in our formulation. 
This helps us to connect our formulation with the previ- 
ous ones [24, 25, 30]. Note also that the memory kernel 
in Eq. (103) cannot be expressed solely in terms of the 
time difference t — s, even after the application of the 
mode-coupling approximation. This is in contrast to the 
cases of CK [26], where the memory kernels depend only 
on the time difference t — s at the level of the Mori-type 
equations. 



VI. STEADY-STATE PROPERTIES 

In the previous sections, we have derived a set of closed 
equations for the time-correlators 3> q (t) and H q (t), i.e. 
Eqs. (55) and (83), where the memory kernels are given 
by Eqs. (102) and (103). According to the ITT scheme 
[24, 28, 30], the steady-state properties are written in 
terms of the time-correlators, with equilibrium quanti- 
ties (e.g. the static structure factor) as the only inputs. 
We follow this scheme and derive a closed formula for 
the steady-state quantities. In the course of this proce- 
dure, we will derive an explicit expression of the multi- 
plier X a {t) introduced in section IV D. A specific formula 
for the steady-state shear stress is derived in section VI C. 



A. Time-correlators of interest 

From the analogue of the Green-Kubo relation 
Eq. (19), time-correlators of interest are of the follow- 
ing form, 

G AB (t) = (A(T(t))B(T) ) = ( [U(t)A(T)} B(T) ) , (104) 
B(T) = <J xy {T) or a(T)5K(T). (105) 

In this paper, we concentrate on the steady-state kinetic 
temperature and the shear stress; i.e. we consider the 
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specific cases, A(T) = K(T) and a xy (T). 

Similarly to the case of static projection operators [26] , 
we can show that the time-correlator Eq. (104) resides in 
the subspace orthogonal to the density and the current- 
density fluctuations, i.e. 



G A , B (t) = {[U (t,0)A(r)]B(r)), 



(106) 



whose proof is given in Appendix C 8. 

Now we apply the second projection operator to 
Eq. (106). Since B(T) = a xy (T) and a(T)6K(T) are 
zero-wavevector quantities, it is sufficient to project onto 
the "zero-mode" pair-density correlator, 



(107) 



V$(t)X = ^ «*(*)"*(*)' 



fc>0 



J k(t) 



V%(0)U o (t,0)-p°(t)A(T)\ B(T) 
Uo(t,0)V$(t)A(T)] V°(0)B(T) ) , (108) 



which is a restricted form of Eq. (84) . The mode-coupling 
approximation to Eq. (106) is then 



G AtB (t) 



where the Hermiticity of the projection operator V® is 
applied in the last equality. Note that U o (t 7 0), rather 
than Uo(t, 0), appears due to the insertion of V^t). 



B. Isothermal condition in MCT 

To proceed further, it is necessary to calculate the ex- 
plicit form of the multiplier X a (t) given by Eq. (80). For 
this purpose, we specify the form of a(T) and explicitly 
calculate GK,<r{t) and GK,aSK{t), which are given in the 
mode-coupling approximation as 



GkA*)^( Uo(t,0)V°(t)K(T) P 2 °(0K y (r)\, (109) 



G 



K,a5K 



(t) 



U (t, 0)VUt)K(T) V$(0)a(T)6K(T) 



(110) 



In the following, we examine concrete examples of a(T) 
and X a (t). The simplest choice of a(T) is to fully ne- 
glect its Independence. In this case, no fluctuations are 
incorporated, and a(T) is a constant, which we denote 
«o- However, we illustrate here that this choice always 
causes a heating- up of the kinetic temperature [39] . For 
a(T) = a , Eq. (110) reads 



GK,aSK(t) 



U (t, 0)V^(t)K(T) V^(0)6K(T) ) , (111) 



From the above argument, it is shown that, in order to 
retain GK,aSK{t) non- vanishing, it is at least necesssary 
to introduce fluctuations into a(T). The simplest choice 
for this is to incorporate current fluctuations as 



which vanishes due to V® (0)8K(r) = 0. This implies, 
together with the fact that Gx,a(t) ^ which is shown 
below, that the isothermal condition Eq. (77) cannot be 
satisfied. Hence, the average kinetic temperature calcu- 
lated this way is unphysical and cannot be regarded as a 
steady-state temperature. Note that this case is the one 
adopted in CK [26]. 



a(T) 



a 



lNk B T 



Pt_ 

^ 2m' 



(112) 



where a is a constant which gives the initial strength 
of the dissipative coupling to the thermostat. Now 
we calculate the projected properties V® (0)[a(T)6K (T)], 
V$(0)a xy (T), and V$(t)K(T) which appear in Eqs. (109) 
and (110), with the choice of Eq. (112) for a(T). Straight- 
forward calculation leads to 



V°(0)[a(T)5K(T)} = a ^ ^ ^ 



N 



E 



fc>0 



n k(t) n k(t) 



1 £?a 5fe « 



(113) 
(114) 



and 



n{t)a xy {T) 



k B T \ - W fc ( t ) 

~~ N h s m 

k x k y 1 dSk 
k Sk dk ' 



n k(t)n* k(t) , (115) 



(116) 



is derived in Appendix C9. From Eqs. (109), (110), and 
(113)— (115), Gx,a{t) and GK,aSK{t) are given by 



GkAQ 



3 (k B T) 2 
2 N 



EE 



1 w k , 

Sk' 



fc>0 fc'>0 ' k ^ 
Uo{t,0)n k{t) n* k{t) ^ n kl n* k , 



(117) 



and 



^K,a5K\t) — -CXq — 



EE 



i i 



fc>0 fc'> 



fc(>0 S k(t) S k ' 



U (t, 0)n fc (t)<( t )J n k >n%, j , (118) 

respectively. Application of the factorization approxima- 
tion to the four-point function reads 



1 

N2 



Uo(t,0)n k{t) n* k(t) 



n k m k , 



[Sk' ,k + S k > ,-k] ®k(t) 2 



(119) 



where n* k = n- k , <fr-fe(t) = $fe(i) are applied. The 
derivation of Eq. (119) is shown in Appendix C 7. From 
Eqs. (80) and (117)-(119), the functional form of the 
multiplier X a (t) should be given by 



X a (t) 



7 



'k>0 S k(t) S k 



<H(t) 5 



2a y: 



k>0 S k(t) S k 



(120) 



in order to satisfy the constraint Eq. (77). This can be 
regarded as an analogue for MCT of the GIK thermostat, 
Eq. (3), for molecular dynamics. 
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We note here that, with the choice of Eq. (112), the 
second term of A(r) inEq. (9) reads ^^daiT) / dpi]-pi = 
2a (r). This is negligible compared to the first term of 
Eq. (9), which is proportional to N, in the thermody- 
namic limit. 



C. Steady-state stress formula 

We define the steady-state shear stress by the steady- 
state stress tensor as follows, 



(<?) 



ss 



(<Txy) 



ss ■ 



(121) 



From the steady-state formula Eqs. (19) and (20), and re- 
minding that (a xy (T) ) — 0, the steady-state shear stress 
is given by the time-correlators as 



ss 



Pi 

V jo 



dtG a>a (t) 



2/3 
V .,„ 



dtGl X) agK (t),(122) 



where the time-correlators G a ^{t) and G^ a5K (t) are 
given by 



(123) 



and 



G 



<j,a5K 



(t) = ([U R (t)a xy (T)]a(T)SK(T)) 

= K{t)G^ aSK {t). (124) 

In Eq. (124), we have introduced the multiplier X a (t) 
by applying U R (t), Eq. (78), which is required from the 
arguments of section IV D. In the mode-coupling approx- 
imation, G ata {t) and G^ aSK {t) are approximated as 

G a ,„{t) ~ ( [u (t,0)V$(t)a xy (T)\ V%{Q)a xy {T)) (125) 
and 



G 



(A) 

cr.aSK 



(*) 



\ a (t) ([u (t,o)v°(t)<j xy (r) 

xP 2 » [a(T)5K(T)]), 



(126) 



respectively. From Eqs. (113), (115), (119), (125), and 
(126), Eq. (122) reads 



k R T 



ss 



dt 



2 L 

d 3 k W k{t) jW k - 2\ a (t)a 



(2tt)3 S m 



$ fe (t) 2 .(127) 



We can see that a correction due to the dissipative cou- 
pling, which originates in the current fluctuation incor- 
porated in a(T), arises to the well-known formula for 
the overdamped case [30]. This is in contrast to the 
case of a(T) = a , where the steady-state stress for- 
mula is coincident with the overdamped case [26]. Note 
that this correction parallels the "cooling effect" of the 
kinetic temperature, which compensates the "heat-up" 
by the shearing and keeps the kinetic temperature un- 
changed from its initial equilibrium value. 



VII. NUMERICAL CALCULATION 

To demonstrate the validity of our formulation, we 
show the results of the numerical calculations in this 
section. As is well known, it is ineffective at present 
to perform grid calculations for a three-dimensional 
system, due to limitations of computational resources. 
In this work, we adopt the "isotropic approximation" , 
which has been formulated and implemented by FC 
[28, 30] and CK [26]. The grid calculations in two- 
dimensional sheared Brownian systems [40, 41] show that 
the anisotropy is relatively small, which assures the vali- 
dation of this approximation, at least in two dimensions 
[42] . This approximation also enables us to compare for- 
mally similar equations for the underdamped and over- 
damped cases, as is discussed in Appendix B. The details 
of the formulation of the isotropic approximation is well 
described by CK [26] , so we only show the results below 
and give some additional remarks in Appendix D. 

Sheared systems are genuinely anisotropic, which can 
be seen from, e.g. the existence of the anisotropic 
term — [n ■ H q (t)] X in the Mori-type equation, Eq. (83). 
By the application of the isotropic approximation, the 
anisotropic terms are neglected, which allows us to ob- 
tain a single second-order equation for the density time- 
correlator by combining Eqs. (55) and (83). The resulting 
equation (MCT equation) is 



dt 2 



„2 m 2 

\(t) 



\ a (t)a - 7 



it 



i + Uity 



r d 

-J dsM q{s) (t-s) — $,(«), (128) 



where the memory kernel is given by 



nvj 



i + \(it) 2 



J (27T 



3 fc 

) 3 



(<? ' &) Ck(s) + (flf ' P) Cp(s) 



(q ■ k) c- k{t) + (q-p) c p{t) 

*£(»)(*- *)*?(«)(*-*)■ 
(129) 



Here, p = q — k is assumed. The notation and the deriva- 
tion of Eqs. (128) and (129) are shown in Appendix D 1. 
Since the multiplier X a (t) given by Eq. (120) is a func- 
tional of <fr q (t), Eq. (128) is solved by iteration between 
<&q(i) and X a (t). 



A. Time-correlators 

In the isotropic approximation, the MCT equation 
Eq. (128) is numerically solved on a one-dimensional 
spatial and temporal grids. The spatial grid is for the 
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wavenumber (modulus of the wavevector). The dis- 
cretized form of the memory kernel Eq. (129) on the 
spatial grid is shown in Appendix D2. For the time 
integration, we adopt the algorithm of Ref. [43], which 
enables us to calculate robustly in long time-scales by 
the gradual coarse-graining in the temporal grid. 

For the units of non-dimensionalization, we choose 
the diameter d and the mass to of the sphere for the 
length and mass, respectively, and tq = d/vo, where 
vq = Wq + ' — Vq = jL is the relative shear velocity at 
the two boundaries, for the time. 

The conditions of the calculation are as follows. The 
spatial grid is chosen as qd = qA, where A = 0.4 is the 
grid spacing and q = (2m — l)/2 (to = 1,2, ••• ,100) 
is the discretized index. The cut-off of the wavenum- 
ber is g max rf = 39.8. The number of the temporal grids 
is Nt = 256. The time-step is initially Ato = 10 -6 To, 
which is doubled in every Nt/2 steps. There are three 
inputs; the volume fraction ip = 7rnd 3 /6, the static struc- 
ture factor S q , and the shear rate 7. The volume fraction 
is expressed in terms of the "distance" from the criti- 
cal volume fraction of the MCT transition in the equi- 
librium MCT, tp c = 0.51591213 [44], which is denoted 
as e e (ip — ip c ) /(p c . This definition of the "distance" 
implies e > for the glass phase, while e < for the 
liquid phase. The value of e is fixed at e = +10 -3 
for the calculation of the time-correlator, while it is 
varied for the shear stress, as is shown in subsection 
VII B. As for the static structure factor, the analytic 
solution of the Percus-Yevick equation [35] for three- 
dimensional hard-sphere systems is adopted, whose ex- 
plicit expression in the Fourier space can be found in e.g. 
Ref. [45]. The initial conditions are $> q (t = 0) = S q and 



m q {t)/dt] t 



0. 



The result of the calculation is shown in Fig. 1. Here, 
the wavenumber is fixed at qd = 7.0 (the first, highest 
peak of the static struture factor; refer to Fig. 4), while 
the shear rate 7 is varied. The lines correspond to 7T0 = 
(no shear), 10" 8 , 10~ 6 , 10" 4 , 10" 2 , and 1. 

For comparison, we show the result for the overdamped 
case in Fig. 2. The MCT equation for this case is 
shown in Appendix B, Eq. (B6). The effective friction 
coefficient which appears in Eq. (B6), a d, is fixed as 
ot dd 2 1 \v\tq) = 0.1. The result of Fig. 2 is qualitatively 
in accordance with the previous results [40, 41], and al- 
most coincident with them at the quantitative level as 
well. 

From Figs. 1 and 2, the density time-correlator de- 
cays due to shearing around the time-scale r Q ~ 7" 1 
(a-relaxation time), for both the underdamped and the 
overdamped cases. The resemblance between the two 
cases can be seen not only in the a-relaxation time r a , 
but also in the non-ergodic parameter (NEP), which is 
almost coincident. These results are consistent with the 
observation that the long-time dynamics after the early- 
/3-relaxation time rp is dominated by the memory kernel, 
and the instantaneous dynamics are invalid at this time- 
scale [31-33]. The difference between the two cases can 
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FIG. 1. (color online) Numerical solution for the normalized 
density time-correlator. The wavenumber is fixed at qd = 7.0 
(the first, highest peak of the static structure factor), while 
the shear rate 7 is varied. The lines correspond to 770 with 
(a) no shear (zero), (b) 10" 8 , (c) 10" 6 , (d) 10" 4 , (e) 10" 2 , 
and (f) 1. 
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FIG. 2. (color online) Numerical solution for the normalized 
density time-correlator in the overdamped limit. The condi- 
tions and the captions are the same as Fig.l. 



be seen at the early stage before 173 [40]. As for the un- 
derdamped case, the density time-correlator is held con- 
stant at its initial value until t ~ 10 To, since the fre- 
quency of the sound wave, u) q {t) = ^Jv^q(t) 2 / S q (t), dom- 
inates the transient behavior at this stage (for qd = 7.0, 
w 9 (t) _1 ~ 0.27r at t <C rg). On the other hand, for the 
overdamped case, the density time-correlator is already 
decreasing at t ~ 10~ 4 to, which can be seen from its ap- 
proximate solution at this stage, & q (t) ~ exp [— t/r d], 
where r d = a dS q ( t )/ [vrlit) 2 ] = a d/^ q {t) 2 is the 
time-scale of this damping (for qd = 7.0, r Q d — 7x 10~ 3 t 
at t <C Tp). The emergence of two time-scales is one 
of the significant features of the underdamped systems. 
In overdamped systems, there is only a single time-scale 
which is a ratio of a d and w q (t), while X a (t)ao and 
u>q(t) settle independent time scales in the underdamped 
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case. Due to this fact, overdamped systems are scaled 
by a single non-dimensional parameter, the Peclet num- 
ber Pe = 7T0, while this is not the case for the under- 
damped case. There are also effects of the difference on 
the steady-state shear stress, which will be discussed in 
subsection VII B. 
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FIG. 3. (color online) Numerical solution for the normalized 
density time-correlator. The shear rate is fixed at 7T0 = 10 -2 , 
while results for several wavenumbers below the first, highest 
peak of the static structure factor are shown. The lines corre- 
spond to wavenumbers qd with 5.0, 6.2, and 6.6. "u.d." and 
"o.d." in the caption correspond to the underdamped and the 
overdamped cases, respectively. 
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FIG. 4. (color online) The structure factor used as an input in 
the calculation. The three wavenumbers whose density time- 
correlator is shown in Fig. 3 are highlighted in (red) solid 
circles; qd = 5.0, 6.2, and 6.6, from left to right, qd — 7.0 
corresponds to the first, highest peak. 



Next, we show the results for several wavenumbers 
in Fig. 3. The shear rate is fixed at 7T0 = 10~ 2 , and 
other conditions are the same as those in Fig. 1. Three 
wavenumbers below the first, highest peak of the static 
structure factor are chosen; qd = 5.0, 6.2, and 6.6. They 
are depicted in Fig. 4 in (red) solid circles, where the 
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FIG. 5. (color online) Numerical solution for the normalized 
density time-correlator of the theory by Chong and Kim [26]. 
The conditions are the same as Fig. 3. 



static structure factor we adopt is shown. We can see 
that the density time-correlator is almost monotonically 
decreasing, except for the spike around t ~ 10 -1 To, for 
the underdamped case. This spike is the vestige of the 
oscillation of the sound wave, which is smeared out in 
longer time-scales. In fact, there is no spike in the re- 
sult for the overdamped case, which is shown together in 
Fig. 3. 

Next, we show the result for the CK theory [26] in 
Fig. 5, where the conditions are the same as those in 
Fig. 3. The difference with the result of our formu- 
lation is obvious; there are significant signals of over- 
shoot/undershoot, i.e. the normalized density time- 
correlator exceeds 1.0 or becomes negative. For equilib- 
rium systems, it is easy to prove that the absolute value 
of the normalized density time-correlator is less than 1.0 
[35], and that the density time-correlator monotonically 
decays in the overdamped limit [2], On the other hand, 
for general nonequilibrium systems, there seems to be 
no rigorous proof of the bounded property or the mono- 
tonicity of the density time-correlators. However, it is 
natural to expect that these properties also hold as well, 
at least for cases with small shear. In addition, it is 
obvious that the overshoot/undershoot is not the result 
of the oscillating nature of the underdamped system. 
The overshoot/undershoot appears at the a-relaxation 
regime, where the instantaneous oscillation is sufficiently 
damped already. From these considerations, we conclude 
that the overshoot/undershoot found in CK theory [26] is 
an artifact of the inappropriate definition of the density 
time-correlator. The problem of overshoot /undershoot 
will be discussed further in section VIII around Eq. (135). 



B. Shear stress 

Now we present the result for the steady-state shear 
stress in unit of fceT/d 3 , where d is the diameter of the 
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sphere, in Fig. 6, which is calculated from the solution of 
the density time-correlator by Eq. (127). The conditions 
are the same as those in section VII A, aside from two 
exceptions. One is the strength of the thermostat in the 
overdamped case a d, which is fixed a dd 2 / (v^to) = 1 
here. This value is chosen to conform with the previous 
work [24], which is the direct reference of our calcula- 
tion. Another is the volume fraction, where four cases, 
e = ±10~ 2 , ±10~ 3 , are considered for underdamped and 
overdamped cases, respectively. The results for the un- 
derdamped case are shown in solid lines, while those for 
the overdamped case are shown in dotted lines. 

As discussed in section VII A, underdamped systems 
are not scaled by a single parameter, the Peclet number 
Pe, in contrast to overdamped systems. Hence, we choose 
as the horizontal axis the non-dimensionalized shear rate, 
fro- 
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FIG. 6. (color online) Numerical result for the steady-state 
shear stress, in unit of fcaT/d 3 . The solid lines are for the 
underdamped case, while the dotted ones are for the over- 
damped case. The four lines for each case are for e = ±10 -2 , 
±10 -3 , where e = (tp — tp c ) /tp c is the distance of the volume 
fraction tp from the MCT transition point ip c . "u.d." and 
"o.d." in the caption correspond to the underdamped and 
the overdamped cases, respectively. 

The results for the overdamped case reproduces those 
found in Ref. [24]. We can see from Fig. 6 that the un- 
derdamped case shows similar tendencies with the over- 
damped case. That is, in the liquid phase with e < 0, the 
shear stress shows the Newtonian behavior, a xy oc 7, for 
small shear rates. For large shear rates, the above linear- 
ity is broken, which signals the "shear-thinning" . In the 
glass phase with e > 0, the shear stress remains finite in 
the limit 7 — > 0, which is nothing but the yield stress. 

At the quantitative level, however, discrepancies can 
be observed. For high shear rates, e.g. 7T0 > 1.0, the 
shear stress is larger for the underdamped case. This is 
since the density time-correlator is held constant in the 
short-time regime t < O.Itq due to the inertia effect in 
the underdamped case, while it is already decreasing in 
the overdamped case, as previously discussed. 
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FIG. 7. (color online) Numerical result for the accumulated 
steady-state shear stress (a(t)}, Eq. (130), in unit of ksT/d 3 . 
The two lines are for the underdamped and overdamped cases, 
respectively, for e = +1CP 3 and 770 = 10 -4 . The steady-state 
value of the shear stress is attatched to each line. 
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FIG. 8. (color online) Numerical result for the mutiplier 
\ a (t)- We plot X a (t)ao in unit of 7. The conditions are 
the same as those for Fig. 7. 



On the other hand, for low shear rates, the yield stress 
which emerges in the glass phase (e > 0) is smaller for 
the underdamped case. This is due to the analogue of 
the "cooling effect", which is explained below Eq. (127). 
To convince this statement, the accumulated shear stress 
as a function of time, which we define by 



< *(*)> = 



k B T f* 
— h dS 

X J si Ms) ' (130) 

is shown in Fig. 7, and the multiplier X a (t), given by 
Eq. (120), is shown in Fig. 8, for the case e = +10~ 3 and 
7T = 10 -4 , respectively. In Fig. 7, the steady-state value 
of the shear stress is attatched to each result. We can 
see that the discrepancy between the underdamped and 
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the overdamped cases arises in the a-relaxation regime, 
t ~ 7 _1 = 10 4 r . It is notable that the multiplier X a (t) 
correspondingly magnifies in this regime, which suggests 
the growth of the current fluctuation incorporated in 
a(T). This implies that the growth of the current fluctu- 
ation in the dissipative coupling relaxes the shear stress 
in the underdamped case. It might also be worth not- 
ing that the density time-correlator are almost coinci- 
dent in the a-relaxation regime for the underdamped and 
the overdamped cases. However, in the underdamped 
case, the time derivative of the density time-correlator 
is nothing but the density-current cross time-correlator, 
Eq. (B2), which also exhibits a growth in the a-relaxation 
regime. This is another evidence that the current fluc- 
tuation incorporated in a(T) is the origin of the stress 
relaxation. 

Now it is clear why the correction in Eq. (127) ap- 
pears in the vertex function; since it originates in the 
current fluctuation, it cannot appear in the density time- 
correlator, which leaves the vertex function as the only 
possibility. Together with the fact that the correction 
is actually the multiplier Eq. (120) which is written in 
terms of the density time-correlator, Eq. (120) can be 
regarded as a current fluctuation, which expectedly is 
highly nonlinear. 

It is remarkable that the existence of a discrepancy 
between the overdamped MCT and the molecular dy- 
namics (MD) simulation, which has an inertia effect in 
the a-relaxation regime, has been reported [46]. Our re- 
sult presented here suggests that the discrepancy is quite 
generic. Thus, we should be careful in comparing the 
results of MD with overdamped dynamics such as Brow- 
nian dynamics or the overdamped MCT. 



C. Response to a perturbation 

As an application of the underdamped MCT so far con- 
structed, we calculate the response of the density time- 
correlator to a perturbation of the shear rate, and demon- 
strate the significance of the underdamped formulation. 

We discuss a response to an instantaneous change in 
the shear rate 7 at the quasi-steady state which corre- 
sponds to the plateau of the density time-correlator at 
t = t . Specifically, we consider a pulse-like perturbation 
of the form 

7(t) = 7 + A 7 [9(t - to) - 9(i - (t + Aty))] , (131) 

where Atj is the width of a rectangular pulse, 0(t) is 
the Heaviside's step function, and to is of the order of 
the early- /3-relaxation time, to ~ Tp. 

For a time-dependent external field (the shear rate in 
our case), the time-evolution operator is expressed in 
terms of a time-ordered exponential [27], and hence we 
cannot simply apply the MCT formulated in this paper, 
which is valid for a constant shear rate. However, for a 
pulse-like perturbation, it can be shown that the time- 
ordered exponential reduces to a normal exponential for 



the case of a weak and instantaneous perturbation, i.e. 
A7/7 -c 1 and Aty/(t — t ) <C 1. In this case, the time- 
evolution operator U->(to,t) can be expanded as 



EM*o,t) 



JC(Aj)(t-t ) 



+ 0(A 7 2 (t-t o ) 2 ), (132) 



where iC(Aj) is the perturbed Liouvillian, whose shear 
part i£j 
given by 



•yiCj, i£j = J2fLi(yi d / dx i-P V t d / d Pf) is 



iC^(Aj) = j{t)i£j, 
7 (t) =7(1 



(133) 



^^-e(t-(fo + A^))|.(m) 

7 t — to 



That is, we can apply the MCT formulated for a normal 
exponential, together with a time-dependent shear rate, 
Eq. (134). 

The conditions of the calculation are as follows. The 
time when the perturbation is switched on, i.e. to, is set 
to the order of the early-/?- relaxation time Tp ~ 10 2 To- 
This choice is made by the observation that the onset 
of the plateau of the density time-correlator is around 
10 2 To, which can be seen from Fig. 1. The shear rate and 
its magnitude of perturbation are set to 7 = 10~ 4 To and 
A7/7 = 0.1, respectively. Note that the a-relaxation 
time is r a ~ 7 _1 = 10 4 t , and hence the quasi-steady 
state lasts for a time interval of 99to- The initial time 
step is set to Ato = 10~ 6 to, which is doubled in every 
4096 steps (N t = 8192). The width of the pulse is set to 
Atj = At , which is small enough compared to the time- 
scale of the response. Hence, the shear rate Eq. (134) is 
pulse-like at this time scale. 

The result for the response of the normalized density 



time-correlator, — 



(A 7 ), 



(t) - 5>q 0) (t) /S q , is shown in 



Fig. 9, both for the underdamped and overdamped cases. 



Here, <&g A ^(t) and &q U '(t) are the perturbed and the un- 
perturbed density time-correlators, respectively, and the 
sign convention is chosen simply for convenience. We can 
see that, in the underdamped case, there is a delay in the 
emergence of the response, and an oscillation of period 
typically of the order of 5 x 10~ 3 to can be seen. This 
feature is due to the inertia effect of the underdamped 
system. On the other hand, in the overdamped case, 
the response emerges instantaneously and then decays 
with the relaxation time 10 _3 to- Hence, even though the 
density plateau coincides for the underdamped and over- 
damped clear difference can be observed in the 
response at this quasi-steady state. 

The result of this subsection can be further applied to 
the response theory. A formulation and calculation of 
the response of the shear stress to a perturbation of the 
shear rate is presented in Appendix A. 



VIII. DISCUSSION 

In this section, we first compare our work with the 
previous works. To the best of our knowledge, the major 
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FIG. 9. (color online) Response of the density time- 
correlator. The upper (lower) figure is for the underdamped 
(overdamped) case. The four lines correspond to non- 
dimensionalized wavenumbers qd = 0.6, 7.0, 7.4, 7.8 



representative works in sheared MCTs are Fuchs-Catesl 
(FC1) [24, 28], Fuchs-Cates2 (FC2) [30], Miyazaki et al. 
[25, 40], Chong-Kim (CK) [26], and Hayakawa-Otsuki 
(HO) [21]. Besides HO, which is formulated for inelas- 
tic granular systems, they are for sheared Brownian or 
thermostatted systems. The relations between the above 
theories might be confusing, so we briefly review their 
basic setups, and then discuss about the resulting formu- 
lations. 

The FC theories, FC1 and FC2, are for overdamped 
systems whose dynamics is governed by the Smolu- 
chowski operator. In FC1, the density time-correlator is 

defined as $f C1) (t) = (n q (t)n* (4) (0) \ /N, with q(t) = 

q + q ■ nt. It is discussed in Ref. [30] that the ap- 
plication of MCT to the above time-correlator leads to 
non-positive-definite "initial decay rates", which causes 
numerical instabilities. This has motivated the modifi- 
cation which leads to FC2. In FC2, the density time- 
correlator is defined as Q q FC2 \t) = ^n q ^ t )(t)n*(0)) /N, 
with q(t) = q — q ■ nt, which is identical to our def- 
inition. Moreover, the resulting MCT equation is also 
formally correspondent if we neglect the transverse mode 



(or adopt the isotropic approximation) and take the over- 
damped limit (i.e. neglect the inertia term) in our frame- 
work; refer to Appendix B for this issue. Hence, our 
formulation can be regarded as an extension of that of 
FC2 to underdamped systems. To be more concrete, the 
initial decay rate T q (t) = q(t) 2 v^/ S q / t ) is properly ad- 
vected, and the memory kernel possesses the structure 
of wavevector indices depicted schematically in Fig. 10; 
i.e. the memory kernel consists of vertex functions with 
wavevector indices advected to time s, V q ( s ) y k(s),k( s )> an< ^ 
time t, V^( t ) fc( t ) jP ( t ), respectively, which are bridged by 
a projection-free propagator starting from time s with 
interval t — s, e.g. &k( s )(t — s). These features seem 
physically sensible, which manifests the "alignment of 
the wavevectors" which we adopt as a principle in sec- 
tions IIIC and IV C. As discussed in Ref. [30], note that 
the application of the time- dependent projection opera- 
tors enables the preservation of this principle. 




time 



FIG. 10. A schematic diagram for the structure of the mem- 
ory kernel in our formulation. The kernels in the theories by 
Fuchs-Cates [30] and Miyazaki et al. [25, 40] also have the 
same structure. 



q(t-s)Mt-s),p(t-s) 



time 



FIG. 11. A schematic diagram for the structure of the memory 
kernel in the theory by Chong-Kim [26]. 

The CK theory starts with an almost identical mi- 
croscopic framework with ours; i.e. an underdamped 
SLLOD equation governed by a Liouvillian. The crucial 
difference at this level is that the dissipative coupling 
to the thermostat is a constant, while it contains current 
fluctuations in order to satisfy the isothermal condition in 
our framework. At the level of the Mori-type equations, 
another crucial difference is that the time-correlators are 




defined as equivalent to FC1, e.g. $^ CK) (i) = ^^'(t) 
for the density time-correlator, which are inequivalent to 
our definition. The conventional static projection opera- 
tors are applied, and the resulting MCT equation differs 
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from ours in two respects: (a) it breaks the alignment of 
the wavevectors, and (b) the memory kernel survives. 
The feature (a) can be seen in (i) the initial decay rate 

Tq 01 ^ = q 2 v^/S q , which is not advected, (ii) the mem- 
ory kernel, which has the structure of wavevector indices 
schematically depicted in Fig.ll, and probably the most 
significantly in (iii) the wavevector-structure of the time 
integration, which reads 

J dsM g (t-s)^ q{t _ s) (s) (135) 

in the isotropic approximation. The advection of the 
wavenumber index appears in the time- derivative of the 
density time-correlator, while it appears in the memory 
kernel in our framework, cf. Eq. (128). Numerically, 
the overshoot /undershoot found in the CK theory is a 
consequence of the term d^^ t ^ s - ) (s)/ds, which shows a 
singular behavior when the advected wavenumber passes 
the first peak of the static structure factor. On the other 
hand, since the memory kernel includes only the density 
time-correlator, not its time-derivative, a singular behav- 
ior is not found in our result. As for the feature (b), the 
memory kernel has been confirmed to be numerically 
negligible [26], at least for the situation in concern, which 
is compatible with our result, = 0. 

In Miyazaki et al. [25, 40], an alternative approach 
is adopted. They start with a generalized fluctuating 
hydrodynamics for the density and the velocity fields 
with a Gaussian noise. Aside from the drawbacks of this 
approach, which are already discussed (e.g. (i) the as- 
sumption of the fluctuation-dissipation theorem which is 
known not to hold in sheared systems, (ii) the formula- 
tion based on the steady-state fluctuations), the result- 
ing MCT parallels that derived by the projection oper- 
ator formalism. The time-correlators are defined in ac- 
cordance with our work and FC2, and it turns out that 
the MCT equation for the density time-correlator in the 
overdamped limit coincides with that of FC2. 

Next, we discuss the novelties which result from our 
framework. First, the relaxation of the shear stress at 
the a-relaxation regime is a genuine feature of our under- 
damped formulation, which cannot be derived from the 
overdamped framework. This is since the above relax- 
ation is due to the growth of current fluctuations, which 
is not considered in the overdamped case. In the over- 
damped case, there might be a possibility to incorporate 
density fluctuations into the dissipative coupling a d, but 
this does not lead to a relaxation of the shear stress dis- 
cussed here. 

Second, as demonstrated in section VII C, we consider 
an application to the calculation of a response. In addi- 
tion, we show in Appendix A the corresponding response 
of the stress tensor, which is formulated by the response 
of the density time-correlator. From these results, it can 
be seen that although there is no difference in the plateau 
of the density time-correlator for the underdamped and 
overdamped cases in the absense of a perturbation, there 
is a clear difference in a short-time response to a pertur- 



bation for the two cases. This is another clear evidence 
for the relevancy to use the underdamped model even 
in glassy systems, where the dynamics is expected to be 
slow and the effects of current fluctuations nor those of 
inertia are believed to be irrelevant. 

Another virtue of the present framework is that its ex- 
tension to other systems, such as an assembly of granu- 
lar particles, is relatively simple. One just has to replace 
the coupling to the thermostat (the term — a(T)pi(t) in 
Eq. (2)) with e.g. the interparticle viscous interactions 
of the Stokes' type. In fact, it is difficult to incorpo- 
rate viscous interactions inherent in granular materials 
into the formulation of FC [24, 28, 30] nor Miyazaki et 
al. [25, 40]. Dissipation has been introduced in MCT by 
HO [21] and Kranz et al. (KSZ) [22], but they seem not 
to be satisfactory. As for HO, which deals with sheared 
granular systems, dissipation is introduced by the inelas- 
tic Boltzmann operator, i.e. the pseudo-Liouvillian [47]. 
The model considered in KSZ is for a driven granular sys- 
tem under a white-noise thermostat, which does not have 
any advected wavenumber. Moreover, the connection be- 
tween the model by KSZ and actual vibrating granular 
systems is unclear. The application of MCT formulated 
here to sheared granular fluids is now in progress [23, 36]. 
We have already found some remarkable features in MCT 
for sheared granular fluids, where the projection onto the 
density-current mode plays an important role in destruc- 
ting the plateau of the density time-correlator [23]. De- 
tails will be reported elsewhere [36]. 



IX. SUMMARY AND CONCLUDING 
REMARKS 

In this paper, we have constructed a nonequilib- 
rium mode-coupling theory for uniformly sheared under- 
damped systems. For such systems, the theory by Chong 
and Kim (CK) [26] has been known, but the discrep- 
ancies with the theory by Fuchs and Cates (FC) [30] 
in the overdamped limit have been left unresolved. We 
have figured out that the formulation of CK is physically 
not sensible; i.e. it does not satisfy the translational 
invariance in the sheared frame, and leads to peculiar 
overshoot /undershoot of the density time-correlator. We 
have performed a reformulation, starting from the re- 
definition of the time-correlators to satisfy the "align- 
ment of the wavevectors" , which is a consequence of 
the translational invariance. The resulting MCT equa- 
tion preserves this alignment, coincides with that of FC 
in the overdamped limit, and leads to a monotonically- 
decreasing density time-correlator. Furthermore, moti- 
vated by the observation that the constant dissipative 
coupling to the thermostat does not lead to a steady-state 
kinetic tcrmperature, we have implemented the isother- 
mal condition at the level of the Mori-type equations. 
It is essential to incorporate current fluctuations in the 
dissipative coupling to satisfy the isothermal condition. 
Hence, it may well be said that we have non-trivially 



19 



extended the FC theory to underdamped systems in a 
physically sensible way. 

Although it has been believed that the equivalence 
of long-time dynamics in underdamped and overdamped 
systems (e.g. the NEP and the scaling properties of the 
density time-correlator at the a-relaxation) also holds 
for sheared systems, we have figured out two deviations 
of the underdamed from the overdamped case. These 
are the relaxation of the shear stress at the a-relaxation 
regime due to a growth of current fluctuations in the dis- 
sipative coupling to the thermostat, and the short-time 
response to a perturbation of the shear rate. These find- 
ings are the genuine features of our underdamped formu- 
lation, which cannot be derived from the conventional 
overdamped framework. We also stress that our for- 
mulation provides a physical explanation to the discrep- 
ancy between MD and the overdamped MCT reported in 
Ref. [46]. 

An attractive feature of our formulation is that it is 
relatively simple to extend to granular systems. In these 
strongly noncquilibrium, nonlinear systems, genuine fea- 
tures of the underdamped dynamics are also expected to 
be observed. The construction of the theory and its nu- 
merical analysis is under investigation, and will be pub- 
lished elsewhere [23, 36]. 
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Appendix A: Application to the response theory 

In this Appendix we apply the MCT to the calculation 
of the response formula for the stress tensor. We discuss a 
response to an instantaneous change in the shear rate 7 at 
the quasi-steady state which corresponds to the plateau 
of the density time-correlator at t = to, where to is of the 
order of the early- /3- relaxation time, Tp. We consider a 
pulse-like perturbation identical to the one specified in 
Eq. (131) in section VII C. For a time-reversible ther- 
mostatted system, 7 _1 is of the order of the a- relaxation 
time, jT a ~ 1. On the other hand, there is a large hierar- 
chy between these two times in the weakly-sheared case, 
Tp <C T a , and hence 7^0 <C 1 holds. In the remainder, 



we perform expansions with 7*0, and neglect terms with 
order 0{j 2 4). 



Formulation 



We first formulate a response formula of the stress 
tensor. We consider a response of the stress tensor at 
t > t + Afcy, 

Sa xy (T(t)) = a^(T(t)), (Al) 

where ai^\T(t)) and a x °y(T(t)) are the perturbed and 
the unperturbed stress tensor, respectively. Our aim is to 
derive a formula for the statistical average of this quan- 
tity at the quasi-steady state. In the absence of pertur- 
bation, the average of a phase-space variable A(T(t)) at 
the quasi-steady state is given by 



(A(T(t))) 



ss 



dT(to)p ss (T,t )A(T(t)) (A2) 



for t > t , where [27, 48] 

Pss(I\to) = p eq (r)exp 



<Mvr(-s)) 



(A3) 



is the distribution function of the quasi-steady state, 
expressed in terms of the initial equilibrium distribu- 
tion function p eq (r) and the unperturbed work function 
Qo(r(t)). The unperturbed work function is given by 

O (r) - -Pi^yHT) - 2f3a(T)6K(T), (A4) 

where we include only the conservative part, 
— PA-yaxy (r), to the the perturbation of the work 
function. While the conservative part of the work 
function is 0(7), the dissipative part, — 2(3a(T)SK(T), 
is 0(^f 2 t). This can be seen from the fact that a(T) is 
accompanied by the time-evolution operator Uu(t) given 
by Eq. (78), where X a (t) is given by Eq. (120), and the 
fact that the isotropic part of k x k v is of the order of 
0(jt). Hence, the following equality holds, 



(A5) 



which validates the approximation mentioned above. We 
should note that jt <C 1 holds, since t ~ to- 

Then, the average of the response of the stress tensor 
at the quasi-steady state is given by 

(Sa xy (r(t))) ss = I dT(to)p^\T,t + At^Scr xy (T(t)) 

+0{A 1 2 ) (A6) 

for t > t + Aiy, where 

P (A ^(T,to + At^^p ss (T,to) 
A 7 ' At 



1- 



l|^jf 7 ds4°J(T(to-s))\ (A7) 
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is the perturbed distribution function at t — to + At^. 
Here, A7 is the strain generated by the perturbation, 



A7 = A^Ai 



7> 



(A8) 



which we fix at a small value A7 -C 1, and take the limit 



Atj later. 



The linear response of the stress tensor is given from 
Eq. (A6) by 



d 



dA 1 



(6a xy (T(t))) 



HS 



A7=0 



A 

At 



J dT(t )p ss (T,t )Sa xy (T(t)) 



/ dsa^(T(t -s)) 
Jo 



AU 



dT{t ) 



dsp ss (T,t - s) 



xSa xy (T(t + S ))a^(T(t )) 



(A9) 



where in the final equality we utilized the fact that the 
pulse width is sufficiently small, At^ <C to, which val- 
idates the following approximations, /?ss(r,io — s) ~ 
pss(T,t ) and 6a xy (T(t + s)) ~ 6a xy (T(t)) (0 < s < 
Aty). Note that Eq. (A9) is equivalent to the conven- 
tional linear response formula [49] if the steady state is 
sufficiently close to equilibrium. 



2. Steady-state distribution function 

The response formula Eq. (A9) is still merely a formal 
expression, and additional approximations are necessary 
to perform a concrete calculation. In this section we de- 
rive an approximate expression for the distribution func- 
tion of the quasi-steady state, Eq. (A3), under a weakly- 
sheared condition. 

It is notable that pss(r, t ) can be rewritten as [48, 50] 



Pss(I\io) = Poq(r) • exp 

+0 {i 3 4) , 



±(<e_>,. + <e + ). y . to 



(A10) 



where 9± is defined as 0± = J ° dtflo(T(±t)), and 
(•••)-yt ' (' ' ' )-f-o arc * ne "conditioned averages". Re- 
fer to Eqs. (3.2),' (3.3) of Rcf. [48] for their definition. By 
expanding the Liouvillian as iC = iCo + O^), there holds 
T(t) = T (t) + 0(7*), where T (t) = e lCot T{0). Then, 
the two conditioned averages which appear in Eq. (A10) 



can be expanded in terms of 7 as 



/°d s (r2 (r(- s ))) 7:0 

JO 



to 



ds(n (T (- s ))), 



+o{i 2 tl), 

rto 



(All) 



/ ° ds (o (r(s))) 7;to = / ° d s <fio(r (*))),. to 

Jo Jo 

+Otft 2 ), (A12) 

where £1 = 0(j) should be reminded. Note that the 
phase-space point which appears in the conditioned av- 
erages are r (0) and T (t ) for (• • • ) 7;0 and (•••) 7:to , 
respectively. As in Ref. [48] , we can show that 



d s (n (r (- s ))) 7;0 



ds<fio(r (*))>, ito .(A13) 



From Eqs. (All)-(A13), the following approximation 
holds, 

f° ds <fi (r(*))>, ito = f ds (n (T(- s ))) 

Jo Jo 

+O(ff ), (AU) 

and hence 

(e_) 7:0 + (e + ) 7:to = 2 f° dsfhhi-s)) 

Jo 



(A15) 



where (£l (r(— s))) = £l (~f(—s)) should be noted. 
From Eqs. (A10) and (A15), we obtain 

/•to 

pss(r,t ) = Peq(r) • exp 



[ ° dsn (j(-s)) 

Jo 



+0 (jh 



5% 



(Ai6) 



where it should be reminded that f*° <isf2o(7(— «)) de- 
pends on the choice of the phase-space point 7. Eq. (A16) 
might seem similar to Eq. (A3), but actually it is not. 



The factor exp 



L" dsflQ('j(—s)) in Eq. (A16) is merely 



a number and factors out from the phase-space integral, 
which is not the case for the factor exp J Q to dstt (T(— s)) 
in Eq. (A3). 

The remaining task is to evaluate the factor 
exp J *° ds£lo(~f(— s)) . One way to perform this is 

to expand f2o(7(~ t)) in terms of its equilibrium en- 
semble average and the correction terms. By intro- 
ducing the Kawasaki transform •y K [51] of 7, and 
from the Kawasaki-transform property (f2o(r(— t))) y . = 
— (flo(T(t))) yK . , f2 (7(— tj) can be approximated as 



1/2 



(<fi (r(t)))**. - (n (r(t)) 2 ) cqJ 

^-(OoTO)) eq , (A17) 
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where the correction term is expected to be negligible 
compared to the equilibrium average. Furthermore, by 
only retaining the conservative contribution to the work 
function as in Eq. (A5), Qo(r) — — /3 7 criy'(r), the expo- 
nent of the factor is approximated as 

J*° dsShM-s)) = PjJJ ds (4°)(r( s ))) e(j 

+0( 7 2 i 2 ). (A18) 

From Eqs. (A16) and (A18), the distribution function for 
the quasi-steady state can be approximated as 



cq 



pss(T,t ) = Pcq (T) • exp /3 7 J ds (a^(T{s)) 

+0 ( 7 2 i 2 ) . (A19) 

Finally, from Eqs. (A9) and (A19), the approximate for- 
mula for the linear response is given by 

d 



<9A 7 



(^xv(T(t))> 



HS 



-/3 exp 



A7=0 
to 



/? 7 y o d S (4 y )(r( s ))) eq 



x (Sa xy (T(t))a^(T(t Q ))) + O ( 7 2 * 2 ) . (A20) 



cq 



3. Mode-Coupling Approximation 

So far we have derived a response formula in terms of 
an equilibrium two-point function. One way to calculate 
this function is to apply the mode-coupling approxima- 
tion, which we perform in this section. 

From the definition of the response, Eq. (Al), the two- 
point function in Eq. (A20) consists of the following two 
terms, 



and 



4^)(r(i))4°)(r(i )) 



cq 



cq 



(A21) 



(A22) 



By inserting the "zero- mode" projection operator V^t) 
onto pair-density modes, the perturbed two-point func- 
tion in Eq. (A21), for instance, can be approximated as 



U Q {t,t )V Q 2 {t)ai^\T{t ))\ V$(to)*$>(T(to))] cc 

(A23) 

Applying the factorization approximation, Eq. (A23) can 
be expressed in terms of the density time-correlator <£> q (t) 



cq 



1 (fcaTeq) 2 

2 V 



d 3 k W fcW Wfc(t„) g(A7) 



(2tt) 3 Sk(t) S k (to) 



^)(*) 2 -(A24) 



where ^fc^j (t) is the density time-correlator in the pres- 
ence of the perturbation. A similar approximation holds 
for the unperturbed two-point function Eq. (A22) as well, 



with $ 



(A-V)/ 



q (t) replaced by its unperturbed counterpart, 
q (t). The resulting expression reads 
d 



(o), 



<9A 7 



(8a xy (T(t))) 



exp 



ss 

J A7=0 



d 3 k W fc(s0 W fe (o) , a 
(2tt)3 5 fc(s0 S fe fc {S > . 



(2tt) 3 5 fc(i ) Sfc^) 



ds / ds' 



;*o) / a (A 7 ) 

So 



(0) 



fe( 



L)W 2 }- 

(A25) 

Equation (A25) is our stress formula for the isothermal 
sheared thermostat system. 
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FIG. 12. Response of the stress tensor calculated from 
Eq. (A25), normalized by their maximum values. 

The result for the response formula Eq. (A25) is shown 
in Fig. 12. In this calculation, the result for the response 
of the density time-correlator, which is shown in Fig. 9, 
is utilized. As already seen in the result of Fig. 9, an 
oscillation can be observed for the underdamped case, 
while a monotonically decreasing behavior is seen for the 
overdamped case. 



Appendix B: Mori- type equations for the 
underdamped and overdamped cases 

In this Appendix, we discuss the relation between 
the Mori-type equations for the underdamped and over- 
damped [30] cases. 

In the underdamped case, the Mori-type equation con- 
sists of two equations, Eqs. (55) and (83). If we decom- 
pose H q (t) into the longitudinal and transverse compo- 
nents with respect to q(t), 

H q (t) = H%(t)+HT(t), (Bl) 
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the longitudinal component can be explicitly given from 
Eq. (55) as 



(B2) 



By differentiating Eq. (55) with time, and combining it 
with Eq. (83), we obtain 

^%{t) = q{t) ■ H q (t) + q(t) ■ j t H q {t) 



2q-K- H q (t) 



- f dsq(t) A L q (t,s)%(s) 
Jo 

ft d 

-J dsM q {t,s)—<S> q {s) 

- f dsq(t) x M^(t,s)H^(s), (B3) 
Jo 

where we have utilized Eq. (B2). Here, the effective fric- 
tion coefficient a q (t) stands for 



a q {t) = a {t) + 2 tl 



(B4) 



and the scalar memory kernel M q (t, s) is defined as 



M q (t,s) =q(t) X M^(t,s) 



q(sY 



(B5) 



On the other hand, the Mori-type equation for the 
overdamped case reads 

d f^d 

-$,(t) = -n2(t)$,(t) - j dsm q (t,s)-<S> q (s), (B6) 

where Vt 2 q {t) = v^q(t) 2 /[a od S q ( t) ] and m q (t,s) = 
M q (t,s)/a d [25, 30, 40]. Here, a d is the effective fric- 
tion coefficient which is related to the bare diffusion co- 
efficient Do by a d = ksT / (mD ). 

Aside from the inertia term, Eq. (B3) is not coinci- 
dent with Eq. (B6) in three aspects: (i) the transverse 
component H q (t) does not decouple, (ii) there exists an 
additional memory kernel L x {t) 1 (iii) the effective friction 
coefficient a q (t) depends on q and t. 

In situations where we can neglect H q (t) and L q (t), 
a concrete example of which is shown below, Eq. (B3) is 
reduced to 



d 



- f dsM q (t,s)^ q (s). (B7) 
Jo " s 

Although Eq. (B7) is formally coincident with Eq. (B6) 
in the overdamped limit (i.e. upon neglection of the in- 
ertia term), an exact correspondence is lacking due to 



the discrepancy between a q (t) and a d- Still, however, 
it is tempting and meaningful to compare the solution of 
Eqs. (B7) and (B6), which is performed in this work. 

One manifestation to neglect H q (t) and L q (t) is to 
resort to the isotropic approximation, which we adopt 
in this work (refer to Appendix D 1 for details) . In this 
approximation, H q (t) is ignored by definition, and it is 
shown in section V that the memory kernel L q (t) ex- 
actly vanishes in the mode-coupling approximation. Fur- 
thermore, the scalar memory kernel defined in Eq. (B5) 
coincides with that for the overdamped case in this ap- 
proximation. 



Appendix C: Details of the derivations 

Large portion of the details of the calculations can be 
found in previous papers, e.g. CK [26] and FC [30]. In 
this Appendix, we show some details which is specific in 
this work. 



1. Steady-state formula 

One of the goals of this study is to derive a formula for 
the statistical ensemble average of a phase-space variable 
at the nonequilibrium steady state. For this purpose, we 
derive an analogue of the Green-Kubo formula which is 
valid for sheared thermostatted systems. 

To begin with, we adopt the "Heisenberg picture" for 
the statistical ensemble average of a phase-space variable 
A(T), defined as 

(A(T(t))) = J dT Pini (r)A(T(t)). (CI) 

Here, A(T) is time-evolved while the distribution func- 
tion remains at its initial value, p- m i(T) [27]. Since the 
system is at equilibrium with temperature T at t = 0, 
Pini(r) is the Maxwell-Boltzmann distribution 

e -PH (r) 

Pini ^ = Jdre-PHoir)' ( C2 ) 

where H (T) = K(T) + U(T) is the Hamiltonian, K(T) = 
J2iPi/2m is the total kinetic energy, and j3 = 1/ (fcsT). 

The starting point to derive the analogue of the Green- 
Kubo formula is the integral expression for the nonequi- 
librium distribution function p(T, t), 



i(r)-07 / dse 
Jo 

2/3 [ dse- lcU 
Jo 



-iCs 



[p ini (T)a xy (T)} 
[p ini (r)a(T)6K(T)} . (C3) 



Here. 



SK(T) ee K(T) - -Nk B T 



(C4) 



(C5) 
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are the zero-wavevector limit of the shear stress and the 
fluctuation of the kinetic energy, respectively, which to- 
gether constitute the "work function" , 

Q(T) = -(3ia xy (T) - 2(3a{T)5K{T). (C6) 

As explained in section IV D, note that the parameter 
a(r) depends on T, whose specific choice in our work is 
shown in Eq. (112). From Eq. (C3), we obtain the follow- 
ing expression for the nonequilibrium ensemble average 
for a phase-space variable A(T(t)), 

(A(T(t))) = J dT Pini (T)A(T(t)) = J dTp(T,t)A(T(0)) 

= (A(T(0)) )-Pif ds(A(T(t))cT xy (T(0)) ) 
Jo 



-2/3/ ds(A(T(t))a(T)6K(T(0))), (C7) 
Jo 

where the adjoint relation Eq. (15) is utilized, and the 
integrations are assumed to converge uniformly. Differ- 
entiation of Eq. (C7) with time and the assumption of 
"mixing" results in the existence of a steady state in the 
limit t — V oo, which leads to the following formula for the 
steady-state ensemble average: 

(A) ss = (A(T))-(3j rd s (A(T(t))a xy (T)) 
Jo 

roc 

-2(3 ds(A(T{t))a(T)SK(T)) . (C8) 
Jo 

Here, abbreviated notation A(T) = A(T(0)) is adpoted. 
This completes the derivation of Eq. (19). 

2. Fourier transform in the sheared frame 

The derivation of Eqs. (21), (22), and (24)-(26) is 
shown here. Since the phase-space variables are defined 
in terms of the phase-space coordinates, F = {ri,Pi}^ =1 , 
we need to define "field variables" to consider the prop- 
erties with respect to spatial transformations, e.g. trans- 
lational invariance. A field variable for a phase-space 
variable A(T) in the experimental frame is introduced as 

A(r,t) = J2MT(t))8(r- ri (t)), (C9) 

i 

where Ai(T(t)) is a coefficient which depends only on the 
phase-space variables. Here, {r} is a coordinate fixed in 
space. Its equation of motion is 

SLA(r,t)=J2§- t MT(t))Kr-r i (t)) 

i 

= iCAi(T(t))5(r - n(t)) = i£A(r, t) 

i 

= [iCo + iC a + i£j r + iC % ] A(r, t), (CIO) 



where the Liouvillians iCo, iC a , iCj r and iCj p are given 
in Eqs. (11), (13), (27) and (26), respectively. Simple 
observation leads to 



i£j r A(r,t) = -r-K 1 



d 



From Eqs. (CIO) and (Cll), we obtain 



A(r,t). 



(Cll) 



d T d 

dt dr 



A{r,t) =iCA(r,t), (C12) 



where 



iC = iCo + iC a + iC 



7p 



(C13) 



is the Liouvillian which is obtained by subtracting i£j r 
from iC. 

Next we move to the sheared frame (r,t), which is 
defined by the coordinate transformation, 



f = r — (k ■ r)t = [1 — Kt] ■ r 



t = t. 



x - {jt) y 

V 

z 



(C14) 



This is a comoving frame with the stretching of the wave- 
lengths due to shearing. The value of a field variable is 
unaltered by the transformation Eq. (C14), so 



A{r,t) = A(f,t). 
Simple calculation leads to 



(C15) 



d 2 



d_ 
dr 



A(r,t) = -~A(f,i), (C16) 



which implies, together with Eq. (C12), 



i) 



dt 



--A{f,t) = iCA(f,t). 



(C17) 



Hence, the time-evolution generator in the sheared frame 
is i£, Eq. (C13). 

Now we move on to the Fourier space. Fourier trans- 
form in the experimental frame is 



A g (t) 



d d rA{r,ty q r ^^2A l (T(t))e lq ^ 



(t) 



(C18) 



where q is some wavevector. It can easily be seen that 
A q (t) satisfies Eq. (CIO) as well, 



^A q (t) = iCA q {t). 



(C19) 



Fourier transform in the sheared frame is 
A q (t) = J d 3 rA(r,ty q -* 
( dr 



- Jd 3 r det (Jpj A r. '/ 1 '• ' 

/ 



d d rA(r,t)e zq - 



iq(l — Kt)r 



(C20) 



24 



where Eqs. (C14), (C15) are applied. Eq. (C15) holds in 
the Fourier frame as well, so comparing Eqs. (CI 8) and 
(C20), we obtain a relation between the wavevectors in 
both frames, 



q = q (1 Kt). 



Inverting Eq. (C21), 



q = q-(l + Kt) = q(-t), 
q(t) = q — q ■ Kt. 



(C21) 

(C22) 
(C23) 



We verify that the wavevector in the sheared frame is 
Amnc-dcformed by the shear-rate tensor k. We choose 
the signature convention which is referred to as the "for- 
ward advection" in FC [30]. It can also easily be seen 
that A q /_ t \(t) satisfies Eq. (C17), 



iCA q( _ t) {t). 



(C24) 



3. The Mori-type equations 

The derivation of Eqs. (67)-(72), i.e. the Mori- type 
equation without isothermal condition, is shown here. 
Let us start with the basic properties, i.e. the action of 
the Liouvillians on the density and the current-density 
fluctuations, n q and j q . The following equalities can be 
verified by straightforward manipulations. First, for the 
density fluctuation: 

iCn q = [i£o + iC a + i£ % ] E e tq r * , (C25) 



\- . Pi 

= — f 



(C26) 



iC a n q = -a(T) Y,P^QfJ2 ^ = °- ( C27 ) 



iC % n q = -E^' kT -^-E eiq ' Tl = °- ( C28 ) 

From Eqs. (C26)-(C28), the action of it, is obtained as 
itn q = iq ■ j q . (C29) 
Next, for the current-density fluctuation: 



= [i£o + iC a + iC % ] E ^e*« , (C30) 



^oj q A = E(-- — + F >- 

i x 



m 9r, 



_d_ 

dp, 



E- 



2_ e iq-Tj 



E vi 



;„n Pj Pi _|_ \ e iq ri 



(C31) 



id* = -a(T) EPi • £ E ^ e "" J = 

j 



A 

7 P Jq 



8 TT^P 



(C32) 



2> 



dpi ^— J m 



^l(K. pi )V^ = -[K-j q ] A . 



(C33) 

Now we calculate the "correlated part", which is the 
first term on the r.h.s. of Eq. (66). It is projected by the 
rescaled projection operator Eq. (61) as 



U(t)P t e iC lHCj^ t) =U(t)Y^ 



e lC ^£j X qW 



•LI* 



Nu 



E ^ 



fc(t) 

) iV^ nfe « W + ( [^(*)] #W ) /V^« W f ' 



(*)/ 7V V 2- 



(C34) 



The two-point functions which appear in Eq. (C34) are obtained by explicit manipulations. 
First, 



»A?' A 



Q(t) 



Kit))- Wq(t) 



iCn 



k(t) 



n k{t) 



n 



= {.fqit) [ik(t)-j k{t) ]*) = ik(tr ( j A (t) j£ ( * t) ) = 

where Eq. (C29) and the fact that terms with odd number of momentum variables vanish, are applied. 



(C35) 
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Next, 



Jfe(t) 



Jfe(t) 



(C36) 



where the term with i£ vanishes due to odd number of momentum variables. In order to calculate the first term of 
Eq. (C36), it is necessary to specify the explicit form of a(T). As explained in section VIB, we adopt the form of 
Eq. (112), which leads to 

< [^«] #w > ^ - Nv t {«o (i + g|) ^ + « A "} 

in the thermodynamic limit. From Eqs. (C34), (C35), and (C37), the result for the correlated part is 



(C37) 



U(t)P t e iC *HCj$ (t) = iv 2 T ^>-n q{t) {t) a # (t) (t) [k ■ j q{t) (t)]' 



(C38) 



This leads to the first three terms in the r.h.s. of Eq. (67). 

Next, we calculate the "random part". The term without time integration is defined as the "random force" in 
Eqs. (68), (69). This is nothing but the fourth term in the r.h.s. of Eq. (67). The part which requires some 
calculation is the one with memory kernels. Projection with the rescaled projection operator Eq. (61) is 



d 



dt U ^ 



(mem) 



d S U(s)V s e tC l s iCe' tC ^ s U (t,s)e lC ^ t R x 



J ° k l 



NS k 



-n k + 



e lC ^ s i£U (t,s)R* 



■q(t) 



'<?(*) 



K{t) T 



where we introduced the abbreviated notation, 

U (t,s) = e- i ^ a U (t,s)e i/ ^r t . 
Substitution of Eq. (C39) into Eq. (56) yields 



i 

N 



(mem) 



iCU {t,s)Rl 



-fA 

= -[ dsL x q {t,s)$> q {s)- [ dsM^(t,s)H^(s), 
Jo Jo 



^)/^^^) + {[^Mt,s)R x q{t) 



(C39) 
(C40) 



(C41) 



where the memory kernels L^(t,s) and M^(t,s) are the ones defined in Eqs. (70) and (71). These are the last two 
terms in the r.h.s. of Eq. (67). 



4. Insertion of the projection operator Q 



The derivation of Eq. (88) is shown here. Let X be an arbitrary phase-space variable, and consider the following 
expression, 

Ul(tJ)e iC ** X* = Ulitjy^t' [P(t') + Q{t')]X*, 

(C42) 



where Ul(t,t') is the adjoint of Uo(t,t'), 

U^(t,t') = exp. 



■J dse lC l s i&Q(s)e- lC l s 



(C43) 
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Here, the adjoint Liouvillians iC) and iC^ r are defined in Eqs. (29) and (33), respectively. The correlated part of 
Eq. (C42) is 



= E 



(X*n k(n ) 



NS 



k(t') 



(X*n Ht , } ) f 



»*(*') + 



(*') 



(C44) 



The action of Uq (t, t') on £ = n and j is 



exp 4 



~ exp 4 



J dse iC l s i&Q(s)e- iC l 
- J* dse lC l s i&Q(s) [1 + £(«)] e"^ 



= 0, 



(C45) 



where we neglected E(s) in the last step, as explained below Eq. (86) in section V. 

I 

From Eqs. (C44) and (C45), the correlated part vanishes, Similar manipulation for Eq. (71) leads to 
and hence 



or, taking the adjoint, 

e- iC l t 'u Q {t,t') = Q{t')e 
This is the desired equality. 



<<'[/ (M') 



(C46) 
(C47) 



iCtj Q {t,s)R x q{ 



</(*) 



Q{s)U' {t,s)Q{t)R 



■i{t) 



AR 



(C51) 



5. The memory kernels 

The derivation of Eqs. (91)-(94) is shown here. We 
start with the ensemble average of Eq. (70). From the 
adjoint relation Eq. (28), 



iCU (t,s)R x 



■«(*) 



'«(») 



U a (t,s)R x (t) (iCn q(s) ^j ^ 
U (t,s)R x {t) ]n q{s) n). (C48) 
Substituting Eq. (90) into the two terms of Eq. (C48), 
X*^ 

= ( [Q{8)Ufa,8)Q(t)R$ w ] Q(s)X*), (C49) 



where AR q ^ is the modified random force defined in 
Eq. (94). 



6. The vertex functions 

The derivation of Eqs. (97)-(100) is shown here. We 
start with Eq. (97). From the definition of Eq. (69), 



R q(t) n *k(t) n *p(t) 



Q{t)iCj x {t) 



n k(t) n p(t) 



n k(t) n p(t) 



V(t)iCj x 



1(t) 



Uo(t,s)R x {t) 



n k(t) n P (t) 

(C52) 



The first term of Eq. (C52) is evaluated, by the use of 
the adjoint relation Eq. (28) and Eq. (C29), as 



where X = i£,n q t s -\,n q t s -\£l. Here, the idempotency and 
the Hermiticity of Q(s), and the abbreviated notation 
of Eq. (93) is used. From Eq. (C29), the first term of 
Eq. (C48) is projected out by Q(s), which leaves 



- ~ \ h(t 



i£n k ( t ) 

■A -n* * 



l p(t) 



Jq(t) n k(t) 



iCn 



p(t) 



Q(s)U^t,s)Q(t)R q(t) Q(s) 



. (C50) 



(C53) 
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The three-point function in Eq. (C53) can be calculated 
explicitly as 



Next we deal with Eq. (98). From the definition of 



m m 



= 5 x »S q - k . p Nv 2 S p . 
From Eqs. (C53) and (C54), 



(C54) 



q(t) 



n l(t) n *p(t) 



= iS q , k+p Nv 2 [k{t) x S p(t) + P (t) x S k{t) ] . (C55) 

An explicit manipulation of the projection operator in 
the second term of Eq. (C52) leads to 



x 

«/(*) 



n t(t) n *p(t) 



q' 
+ < 



= E{(F^m 



nq ' {t) I NS q , (t) \ n 9'(*) n fcW"p(t) 



Vv ^ \ 'q (t) n *k(t) n *p(t) 



NS, 



Jq(t) 



iCn„i 



<?'(*) 



V(«) n fe(t) n p(t) 



iq'(ty 



WW) — 



(t) / \ n q'(t) n *k(t) n * P (t) 



= iv T 



S q {t) 



l q(t) n k(t) n p(t) / • 



The convolution approximation [35], 

(n q n* k n* p ) ~ S q ^ k+p NS q S k S p , 



(C56) 



(C57) 



is applied for the three-point function, which finally leads 
to the following expression, 

-ftg(t) n fc(*) n p(t) 

= i^fc+pWr [fc(0 A ^ P (t) +p(*) A ^few] 

-fer | ' Sq.k+pN S q (t)S k (t)Sp(t) 

= i5 q , k+p Nv 2 T [k{t) x S p{t) + p(t) x S k(t) - q(t) x S k{t) S p(t) ] 
= -iS q .k+pNv^S k ^t)Sp(t) 



iS q , k +pNv^S k{t) Sp {t) n [k(t) x c k{t) + p(t) x c p{t) ] . 

(C58) 

Here, c g is the direct correlation function [35] which is 

;t( 

1 



related to the static structure factor S q as 



nCq = 1 S 



(C59) 



by the Ornstein-Zernike relation [35]. From Eq. (C58), 
it is straightforward to verify Eqs. (97) and (99). 



n Ht) n pW AR x * {t) 
N 2 Sk(t)S p (t) 



n k{t)n p (t)R q l t) 



k(t)n p(t) Q(t) j x (t) n 



N 2 S, 



N 2 S k (t)Sp( t ) 



(C60) 



where the first term of Eq. (CGO) is the complex conjugate 
of Eq. (97). Hence, we only need to handle the second 
term. There are two terms in the modified work function, 
Q(r) = ~Pja^ n) (r)~2Pa(r)5K(r). For the first term, 



C(*) 



•A ,_(kin) 
3q(t) a xy 



j x qit) ^ n) -nt) 



i X rr( kin ) 
Jq{t) u xy 



(C61) 



where aiy 1 "^ = J2iPfPi /% m is the kinetic part of the 
shear stress. The projected term is 



v{t) [j q A (t) < in) ] = E 



(kin) 

xy 



+E 



NS m 

La (kin) 
[ J q(t) ax y 



3h(t) 



n k (t) 



3 k (t) 



= mt4(«%+*%). (C62) 

is an odd 



■A (kin) 

J q (t)" x y 



From Eqs. (C61) and (C62), Q(t) 

function of the momentum variables, and hence vanishes. 
For the second term, 



G(t) 



A MT)SK(T) 



3q(tV 

= & t) a(T)6K(T)-V(t) j x {t) a(r),)K(T) .(Cm) 



where SK(T) = J2iPi/2m-3Nk B T '/2 is the fluctuation 
of the kinetic energy, and a(T) is given in Eq. (112). The 
projected term is 



V{t) [j x {t) a(T)5K(T] 



E 



[j x {t) a(T)6K(T) 



l fe(t) 



NS, 



k(t) 



E 



3k(t) 



Nv 2 , 



n k (t) 



•?fc(t) 



(C64) 



;A MT)5K(T) 



3q(t) 1 



is an 



From Eqs. (C63) and (C64), Q(t) 
odd function of the momentum variables, and hence van- 
ishes as well. This completes the derivation of Eq. (98). 
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Finally we show Eq. (100). The part which corresponds 
to the first term of the modified work function fi(r) is 



G(*) 



„ ^(kin) 



n q(t) CT iy 



(C65) 



where the projected term vanishes, 



(kin) 
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(kin) 
n q(t) <J xy 



n>k(t) 



Nv 2 , 



3 Ht) I .[i 
3k(t) 



(C66) 



since (pfpf ) p — 0. Here, (--^p is the ensemble aver- 
age with only the momentum variables integrated. From 
Eqs. (C65) and (C66), 



nk{t)n p ( t )Q{t) 



77* , ^ (kin) 

n q{t) a xy 



= \n k(t) n p{t) n q{t) a^ 



(C67) 



which again vanishes due to {p\p\) p = 0. The second 
term is 

Q(t) [n q(t) a(T)SK(T)] 

= n q{t) a(T)SK(T) - V{t) [n q{t) a(T)SK(T)] . (C68) 

Here, the projected term is 

V{t) [n q{t) a(T)SK(T)} 



% [n q(t) a(T)SK(T)] n* (t) 
E NS~- ~ nk{t) 



k(t) 
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'Jk(t) 



a k B Tn q ( t ), 



(C69) 



where 



[a(T)SK(T)} p = a k B T 



is utilized. From Eqs. (C68) and (C69), 



(C70) 



n k (t)n P (t)Q{t) [n q{t ) a ( T ) SK ( T ) 

= (n k{t) n p{t) n q{t) a(T)5K(T) 

~a k B T ( r7 fe(t) 77 p ( t) n* (t) 



(C71) 



which again vanishes due to Eq. (C70). 

7. Factorization approximations 

The derivation of Eqs. (101) and (119) is shown here. 
From the definition Eq. (93), the left-hand side (l.h.s.) of 
Eq. (101) can be expanded as 



N 2 



Uo(t, s)n fe ( t )n p ( t ) 



n k'(s) n p'(s) 



N 2 

1 



e iC ''r s Uo{t,s)e tl -^ t n k(t) n p(t) 



n fc'(s) n p'(s) 



= 772 ([Ua(t,s)n k n p } n* k ,n* pl ) , 



(C72) 



where the property e lCA,rt n k ^n p ( t ) = n k n p is utilized. 
It is clear that this is the propagator for the pair-density 
fluctuations with respect to the projected time-evolution 
operator Uo(t,s). As is familiar in conventional MCT, 
Eq. (C72) is approximated by factorizing into a product 
of propagators for the density fluctuation with respect to 
the projection-free time-evolution operator e 1 ^* -5 ), 



■^2 ( t* 7 ^*' S ) n k n p] n* k ' n p' ) - ^2 



N 2 
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N 2 
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J£\ r s e iC(t-s) e -i£y r t 



n k 



e i£(t-«) e -«^ P (t-») nfcW 



e iC\ r s e iC(t-s) e -iC^ r t 



V 
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e iC(t-s) e -iC^{t-s) n ^ 



V( S ) 



= N2\ - s ) n Hs)] K'(s) ){[U(t- *)«p(a)] 

= *fe',feV,P*fe(«)(* - s )*p(s)(* - s )- 



(C73) 



In the first step of the r.h.s. of Eq. (C73), the advection 



iC s i 

generators e ^ and e 



arc inserted to account for 



the time when the propagation starts, s, and ends, t. 
Hence Eq. (101) is shown. 



Similarly, the l.h.s. of Eq. (119) is 
1 



N 2 



Uo(t,0)n k{t) n* k{t) 
1 

~ N 2 



n k <n k , 
[U (t,0)n k n k ] n k m%,), 



(C74) 
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which is approximated as 
([U (t,0)n k n* k }n k ,n* k ,) 



Only the potential part of the shear stress cxiy *-* survives 
in the correlator of Eq. (C80), since the kinetic part van- 



ishcs due to (pfpl-) = 0. The remaining part is 



1 



g*-^ g — i£"y*- 1 



e" £ *e 



^ v ot) n* k(t) n k(t) ) = Y.( x ^ el)Vn Ht)<(t) ) ■ (C81) 



+^([e iCt n k{t) ] W>4 



TV 



'fe(t) 



= 5fe',fe$fe(t) • Sk',k^-k(t) + 5- k ',k®k(t) ■ 5k>-k®-k{i) 

= [<**',* + 5*',-fc]**(t) 2 . (C75) 

Hence Eq. (119) is shown. Note that Eq. (C75) is a spe- 
cial case of Eq. (C73) with s = 0, aside from the possible 
pairings of the density fluctuations. 

8. The projected time-correlator 

The derivation of Eq. (106) is shown here. From 
Eq. (64), there holds 

U(t)A(T) = U (t,0)A(T) 



+ I dsU(s)V s e c ^ s iCe- lC ^ s Uo(t,s)A(r). 
Jo 

(C76) 

For the second term, the application of the rescaled pro- 
jection operator Eq. (61) leads to the form 



/' 

Jo 



C, 



(n) 



C, 



(i)A 



ds £l^„ k ,, w+ ^ w(s) j, (c77) 

where (£ = n,j) is a correlator whose detailed ex- 
pression is not important for our purpose. The time- 
correlators of the density and the current-density fluctu- 
ations with zero-wavevector variables B(T) vanish, 

(n h(a) {a)B(T))=5 hfi (n k= o(a)B(T))=0, 
[ti w {8)B{T) ) = 6 k , U = o(s)B(T) ) = 0, (C78) 
since n k=0 — holds by definition, and 

jLo(s)=e^J2^=° ( C79 ) 

i 

from the definition of the peculiar momentum, pi = 
m(fi — k ■ ri) . This completes the derivation. 

9. The projected shear stress 

The derivation of Eq. (115) is shown here. The shear 
stress is projected by the second projection operator as 



V° 2 (t)a xy = 



k>0 



T xy n k(t) n Ht) 



nk(t)K(t)- ( C8 °) 



We utilize the relation well known in equilibrium statis- 
tical mechanics [35], 

which also holds in our formulation since the ensemble av- 
erage ( • • • ) is defined by an averaging with the Maxwell- 
Boltzmann distribution, Eq. (C2). Then, it is straight- 
forward to show 



4T )n kt) n Ht) 



n Ht)n k{t) 



= -k B Tk(t) v J2\ lx ^ lk{t) ' r 



l k(t) 
dnt 



-ik(t)-r? 



n k(t) 



-k B Tk{tf ( a^n fcm + -7^u7^ n Ht) 



= -k B Tk(ty 



\dk(t)x fe w dk(ty 

d 



dk(ty 



n k(t)T1 k{t) 



= -Nk B Tk(t) 
= -Nk B T 



y dS k(t) 



dk(ty 

k(t) x k(t)v as m 



k(t) dk{t) ' 
This, together with Eq. (C80), proves Eq. (115). 



(C83) 



Appendix D: Miscellaneouses of the numerical 
analysis 

1. The isotropic approximation 

The basics of the isotropic approximation and the 
derivation of its resulting equations, Eqs. (128) and (129), 
is shown here. Refer to CK [26] for further details. 

The fundamental idea of the isotropic approximation 
is to reduce to dependence of the three-dimensional 
wavevector to its modulus. This is accomplished for the 
time-correlators by the assumption of 



(Dl) 
(D2) 



In addition, there are polynomials of the wavevectors to 
be handled. They are approximated by their mean values 
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with respect to the solid angles. For instance, 



where q(t) is defined as 



q(t) ■ k(t) = qk (jt) (q x k y + q y k x ) + ( A it) 2 q x k x 
1 



(q-k) 



k x (t)ky(t) = k x [ky ~ (jt)k x ] ~ --{jt)k 2 . 



(D4) 



Note that the anisotropic terms are neglected in the 
above approximations, which are the leading terms in 
the shear-rate for the case jb -C 1. 

The modulus of the advected wavevector is approxi- 
mated as 



(D6) 



(D3) Combining Eqs. (55) and (83), with the application of 
Eqs. (Dl) and (D2), leads to 



q{t) 2 = q 2 - 2(jt)q x q v + (jt) 2 q 2 x ~ q(t), (D5) Eq. (B5), reads 



. q x q y (t) d 

7 qitf dr qy 

/"* d 
-J dsM q (t,s)—$ q (s). (D7) 

The scalar memory kernel M q (t,s), which is defined in 



M q (t,s) 



d 3 k 



2q(s) 2 J (2tt) 3 
x$ fc(s) (t-s)$ p(s) (t-s), 



[(q(t) ■ k{t)) c fe ( t ) + (q(t) ■ p(t)) Cp (t) ] [(q(s) • fc(s)) c fe(s) + (q(s) • p(s)) c p{s) ] 



(D8) 



where Eqs. (103) and (99) are utilized. It is convenient to shift the integration variable in Eq. (D8) as k — > k' = k(s), 
which leads to 



M q (t,s) = M q{s) (t-s), 



where 



M<,(t) = 5 ,/ (i^ [(<? ' k) Ck + iq ' P) Cp] [ (<7(t) ' fe(T)) Ck(T) + iqiT) ' P(T)) Cp(r) ] ^k(r)%(r). (D10) 
Application of Eqs. (D3) and (D4) to Eq. (D7) leads to Eq. (128), and to Eq. (D10) leads to Eq. (129). 



d A k 



(D9) 



2. Discretization 

The discrctized form of the memory kernel on the one-dimensional spatial grid is given by 

<-oo r-q+k 

dk dp (q 2 + k 2 - p 2 )c- k(T) + (q 2 -k 2 + p 2 )c p(T) 

JO J\q-k\ l 

X [(q 2 + k 2 - p 2 )c k + (q 2 - k 2 + p 2 )c p ] fcp$ fe (r)$ p (r) 



M q (r) 



1 nv 2 -. 
32tt 2 q 3 



l + \(ir) 2 



nv\ A 5 1 



32tt 2 d 5 q 3 

k V 

X (q 2 + k 2 - f )c k + (q 2 -k 2 + p 2 )c p ] P%(t). 



(Dll) 



In the last step, the wavenumber is discrctized as qd — Aq, where A is the grid spacing, and q is the discretized index 
of the wavenumber, q = (2m — l)/2 (m = 1, 2, ■ • ■ , M), M being the number of grids. The summation with respect 
to p is restricted to those which satisfy the triangle inequality, i.e., 

p=g+fe-l/2 

2> E ■ 

P p=|g-fc|+l/2 
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